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The  purpose  of  this  study  is  to  present  an  efficient 
implementation  for  optimal  structural  shape  synthesis  which 
is  based  on  a combination  of  the  boundary  element  method 
(BEM)  with  nonlinear  programming  (NLP)  optimization 
algorithms.  Shape  optimization  involving  two-dimensional 
elastic,  torsional  shaft  problems,  and  small-strain 
elastostatics  problems,  is  introduced  in  this  study.  An 
approach  for  grid  refinement  and  grid  adaptation  for  the  BEM 
is  developed  for  application  in  plane  elasticity  problems. 

To  prevent  distortion  of  the  original  domain,  master  nodes 
are  introduced  to  control  the  geometry  of  the  structure.  Two 
strategies  for  the  solution  of  the  adaptive  schemes  are 
investigated,  one  based  on  a variational  approach  and  another 
employing  nonlinear  programming  methods.  The  use  of 
nonlinear  programming  methods  in  grid  distribution  allows 
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greater  flexibility  in  selecting  the  form  and  relative 
importance  of  criteria  employed  for  grid  adaptation  and 
refinement.  By  utilizing  the  characteristics  of  both  the  BEM 
and  the  LU-decomposition  methods,  an  efficient  use  of  a semi- 
analytical  sensitivity  analysis  for  optimum  design  is 
obtained. 

This  study  also  examines  the  use  of  BEM  in  optimal  shape 
synthesis  from  the  standpoint  of  an  integrated  formulation. 

A set  of  linear  eguations  obtained  from  the  boundary  element 
analysis  is  treated  as  equality  constraints  in  the 
mathematical  programming  optimization  problem. 

Several  numerical  examples  are  presented  for 
illustrating  the  effectiveness  and  convenience  of  using  BEM 
with  these  techniques. 


vii 


CHAPTER  1 


INTRODUCTION 


1 . 1 Motivation  and  Objectives 

In  most  structural  design  problems,  the  geometric 
configuration  of  the  structure  is  usually  fixed,  and 
structural  member  dimensions,  such  as  cross-sectional  areas 
of  bar  elements,  cross-sectional  properties  of  beam  elements, 
and  thickness  of  membrane  and  plate  elements,  are  changed  to 
obtain  the  optimum  design.  It  is  obvious,  however,  that  an 
additional  mass  reduction  can  be  obtained  by  changing  the 
structural  geometry  in  addition  to  changing  structural  member 
dimensions.  The  problem  of  optimum  shape  design  includes 
design  variables  that  characterize  either  the  shape  of  an 
exterior  boundary  of  the  structure,  or  the  location  and  shape 
of  an  interior  cutout.  An  extensive  survey  of  attempts  in 
shape  optimization  is  presented  in  Haftka  and  Grandhi  [1]. 

There  are  two  major  components  that  must  be  suitably 
integrated  when  attempting  optimum  structural  design.  The 
first  relates  to  choosing  a mathematical  programming 
algorithm  to  achieve  an  optimal  objective  function  and 
simultaneous  satisfaction  of  a set  of  design  constraints. 
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The  second  is  the  selection  of  an  analytical  tool  for 
determining  the  structure  behavioral  response  needed  in  the 
computational  procedure.  In  shape  design  problems  where 
design  variables  control  the  geometry  of  the  structure, 
additional  degrees  of  difficulty  are  introduced.  At  the  very 
outset,  there  is  generally  an  increase  in  computational 
resource  requirement,  as  it  is  more  difficult  to  obtain  good 
sensitivity  information  with  respect  to  shape  design 
variables,  than  with  respect  to  member  sizing  variables. 

This  problem  is  further  compounded  if  the  finite  element 
method  is  used  as  the  analytical  tool.  Since  the  finite 
element  domain  is  changed  continuously  in  the  redesign, 
undesirable  element  distortion  is  often  introduced,  making  it 
difficult  to  ensure  that  the  accuracy  of  the  finite  element 
analysis  remains  adequate  throughout  the  design  processes. 
Although  adaptive  mesh  generation  procedures  such  as  r-  and 
^“finite  element  methods  have  been  developed  to  reduce  the 
error  due  to  element  distortion  [2],  it  is  cumbersome  to 
perform  a new  mesh  generation  for  each  iteration  of  the  shape 
design. 

The  boundary  element  method  (BEM)  has  been  successfully 
employed  in  various  engineering  disciplines  [3],  and  offers 
several  advantages  over  the  finite  element  method  (FEM)  [4], 
Among  these,  the  more  significant  include  a replacement  of 
the  FEM  mesh  generation  for  the  entire  domain  by  a boundary 
discretization  in  the  BEM,  a higher  accuracy  in  sensitivity 
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evaluation  along  the  boundary,  and  simpler  preparation  of 
input  data  for  the  problem.  Furthermore,  the  level  of 
experience  necessary  in  using  BEM  is  generally  not  as  high  as 
is  reguired  in  the  use  of  FEM  for  mesh  generation  and 
selection  of  optimum  element  sizes. 

However,  some  judgement  is  still  reguired  in  using  the 
BEM  for  optimum  shape  design.  In  order  to  reduce  the  number 
of  grid  points  in  the  boundary  element  analysis  and  to  obtain 
a more  accurate  numerical  solution,  a technigue  using  a 
combination  of  grid  refinement  and  grid  adaptation  is 
explored  in  this  study.  In  the  conventional  application  of 
adaptive  technigues,  the  variational  approach  is  usually  used 
to  redistribute  the  analytical  nodes.  In  the  present  study, 
a more  general  nonlinear  mathematical  programming  approach  is 
developed  for  this  purpose.  Such  a strategy  allows 
consideration  of  multiple  criteria  in  the  grid  distribution 
problem. 

The  sensitivity  of  structural  response  to  changes  in  the 
design  variables  provides  valuable  information  for  design, 
and  is  extensively  used  in  most  efficient  nonlinear 
programming  algorithms.  Such  sensitivity  may  be  obtained  by 
either  finite  difference  methods  or  an  explicit  analytical 
differentiation  of  the  response  eguations.  The  former  is 
easier  to  implement  but  has  an  associated  computational  cost 
that  may  be  undesirable  in  large  problems.  The  analytical 
approach  is  sometimes  difficult  to  implement  due  to  the 
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implicit  nature  of  the  response  equations.  In  the  present 
work,  a semi-analytical  approach  based  on  the  characteristics 
of  the  boundary  integral  equation  is  investigated. 

The  concept  of  integrated  approach  in  structural  design 
is  also  studied  in  this  work.  This  approach  was  first 
presented  by  Schmit  and  Fox  [5,6]  in  a FEM  based  structural 
design.  The  basic  idea  behind  the  effort  was  to  treat  the 
behavioral  response  equations  (finite  element  equilibrium 
equations)  as  additional  equality  constraints  in  the 
optimization  problem.  The  behavior  response  variables  were 
added  to  the  list  of  design  variables.  The  equality 
constraints  and  the  inequality  design  constraints  were 
subsequently  handled  by  a penalty  function  approach  in  an 
integrated  design  problem.  The  advantage  of  the  approach 
resided  in  the  fact  that  no  explicit  solution  of  the  response 
equations  was  necessary,  and  the  optimum  solution 
simultaneously  met  the  desired  objective,  satisfaction  of  the 
response  constraints,  and  a solution  to  the  analysis  problem. 
Ineffective  methods  of  solving  equality  constrained  problems 
did  not  allow  this  method  to  be  efficient  till  very  recently 
[7].  In  the  present  work,  a generalized  reduced  gradient 
approach  is  used  for  the  solution  of  the  equality  constrained 
problem  in  such  an  integrated  formulation. 
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1 . 2 Scope  of  Present  Work 

The  literature  regarding  optimum  shape  design  problems 
is  presented  in  Chapter  2.  An  overview  of  mathematical 
programming  methods  and  the  corresponding  mathematical 
formulations  are  also  presented  in  the  same  chapter. 

In  Chapter  3,  the  boundary  integral  equations  of  two- 
dimensional  elasticity  and  torsional  problems  are  derived. 

The  applications  and  numerical  accuracy  of  the  BEM  for  both 
cases  are  illustrated. 

In  Chapter  4,  the  selection  of  shape  design  variables  is 
outlined.  Several  shape  optimization  problems  are  attempted 
by  using  the  BEM  and  a nonlinear  programming  technique  for 
analysis  and  optimization,  respectively.  The  subject  of 
sensitivity  analysis  is  explored  in  Chapter  5.  The  form  of 
the  boundary  integral  equation  is  used  advantageously  in 
deriving  a semi-analytical  approach  for  sensitivity  analysis. 

The  grid  refinement  and  grid  adaptation  techniques  are 
investigated  in  Chapter  6.  Both  the  variational  approach  and 
the  nonlinear  programming  method  are  employed  to  achieve  grid 
adaptation.  Numerical  examples  of  shape  optimization, 
combined  with  these  techniques,  are  also  presented. 

In  Chapter  7,  the  concept  of  the  integrated  approach  is 
extended  to  shape  optimization  problems,  wherein  the  BEM 
analysis  equations  are  represented  as  a set  of  equality 
constraints  in  the  optimization  problems. 
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Finally,  conclusions  drawn  from  this  study  and 
recommendations  for  future  research  are  presented  in  Chapter 
8. 


CHAPTER  2 


OVERVIEW  OF  OPTIMUM  SHAPE  DESIGN 
2 . 1 Literature  Survey 

With  the  increased  availability  of  high-speed  digital 
computers  since  the  1960s,  structural  optimization  has 
occupied  an  important  place  in  the  field  of  engineering 
design.  Recent  studies  in  structural  optimization  have 
extended  the  traditional  structural  member  sizing  problem  to 
the  determination  of  optimal  shapes  for  desired  structural 
response.  The  technique  of  structural  shape  optimization 
has  been  demonstrated  successfully  in  various  publications. 

A general  review  of  papers  published  during  the  last  decade 
can  be  found  in  reference  [l].  In  this  chapter,  the  review 
of  work  related  to  optimum  shape  design  is  classified 
according  to  the  two  analytical  tools  known  as  the  finite 
element  method  and  the  boundary  element  method. 

2.1.1  Finite  Element  Method  for  Optimum  Shape  Design 

The  first  numerical  treatment  of  shape  optimization  was 
proposed  by  Zienkiewicz  and  Campbell  [8]  in  1973.  The 
importance  and  applications  of  shape  optimal  design  was 
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emphasized  in  this  article.  A finite  element  model  was 
formulated  to  analyze  the  shape  optimization  problem,  and 
the  coordinates  of  boundary  nodal  points  were  treated  as 
design  variables.  The  physical  problem  in  this  study  was 
the  design  of  dam  structures  using  isoparametric  elements. 
Sequential  linear  programming  techniques  were  employed  for 
optimization  purposes,  with  the  derivatives  of  structural 
responses  obtained  by  an  implicit  finite  difference  method. 
Ramakrishnan  and  Francavilla  [9]  characterized  this  shape 
optimization  problem  in  more  detail.  A hybrid  penalty 
function  method  was  employed  to  obtain  optimum  design  for 
dam  structures.  The  structural  domain  was  discretized  into 
two-dimensional  isoparametric  elements,  and  the  coefficients 
of  the  element  polynomials  were  chosen  as  design  variables. 
In  1975,  Francavilla,  Ramakrishnan  and  Zienkiewicz  [10] 
first  applied  the  finite  element  method  in  a fillet  shape 
optimization  problem  with  the  objective  of  minimizing  the 
stress  concentrations  in  the  structural  domain. 

Kristensen  and  Madsen  [11]  formulated  a class  of  shape 
optimal  design  problems  for  planar  solids.  They  used 
orthogonal  polynomials  to  locate  the  boundary  of  the 
structural  body  and  treated  the  coefficients  of  these 
polynomials  as  design  variables.  A finite  element  approach 
was  employed  to  obtain  the  structural  response  and  the 
derivatives  of  stress  with  respect  to  design  variables. 

Chung  and  Haug  [12]  used  an  approach  based  on  the 
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calculus  of  variations,  combined  with  adjoint  differential 
equations,  to  obtain  the  sensitivities  necessary  to  perform 
the  two-dimensional  shape  optimal  design.  The  design 
objective  was  to  minimize  structural  weight  with  constraints 
expressed  by  von  Mises  yield  stress  along  the  boundary.  A 
generalized  steepest  descent  algorithm  was  implemented  for 
determination  of  the  optimum  shape.  To  simplify  the 
optimization  procedure,  Dems  and  Mroz  [13]  modified  the 
design  variables  to  a set  of  prescribed  shape  functions  and 
a set  of  shape  parameters;  an  optimality  criterion  was 
derived  and  satisfied  by  using  an  iterative  numerical 
algorithm. 

The  two-dimensional  torque  arm  design  problem  was  first 
solved  by  Botkin  [14],  resulting  in  the  basic  concept  for 
shape  optimization  of  plate  and  shell  structures.  Later, 
Queau  and  Trompette  [15]  used  an  automated  mesh  generator 
for  several  two-dimensional  problems. 

Shape  optimization  of  three-dimensional  structural 
components  was  first  proposed  by  Imam  [16].  In  order  to 
reduce  the  number  of  design  variables  in  the  three- 
dimensional  shape  optimization  problem,  Imam  suggested  the 
use  of  low  order  curves  and  surfaces  to  model  the  domain. 

The  isoparametric  representation  of  the  surfaces  and  the 
numerical  superposition  of  shapes  to  simplify  the  shape 
representation  was  also  proposed.  A feasible  directions 
algorithm  of  nonlinear  programming  was  used  to  obtain  the 
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minimum  mass.  The  dilemma  of  finite  element  distortion  in 
each  new  shape  design  during  the  optimization  procedure  was 
also  discussed  in  this  work. 

The  initial  configuration  of  the  finite  element  mesh  may 
result  in  very  distorted  elements  for  the  new  shape. 
Therefore,  even  the  best  mesh  design  for  the  initial  shape 
is  usually  not  appropriate  for  intermediate  or  final 
designs.  As  a result,  the  analysis  may  be  inaccurate  in  the 
computations  for  stresses  and  displacements.  A solution- 
based  adaptive  mesh  refinement  scheme  was  studied  by  Bennett 
and  Botkin  [17],  in  which  the  strain  density  gradients  were 
used  to  identify  the  regions  which  require  mesh  refinement. 
This  concept  was  also  successfully  extended  to  three- 
dimensional  structural  shape  optimization  [18]. 

Kikuchi  et  al.  [19]  treated  the  automated  mesh 
generation  as  an  integral  part  of  shape  optimization.  Two 
general  grid  adaptive  methods,  the  r-  and  h-  methods, 
applicable  in  shape  optimum  design  are  discussed  in  detail 
in  their  study.  They  emphasize  the  strong  dependence  of  the 
final  shape  of  design  boundaries  on  the  mesh  generation  and 
conclude  that  a good  adaptive  mesh  refinement  strategy  can 
avoid  the  jagged  shape  produced  by  using  the  coordinates  of 
boundary  nodes  as  design  variables.  Furthermore,  an 
accurate  analysis  for  the  discrete  model  may  improve  design 
convergence  and  reduce  the  total  computational  time. 

Using  a mixed  variational  formulation,  Rodrigues  [20] 
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derived  the  optimality  condition  for  which  the  maximum  value 
of  the  von  Mises  equivalent  stress  in  a two-dimensional 
linear  elastic  body  is  minimized.  Choi  and  Seong  [21] 
proposed  domain  integrals  for  shape  design  sensitivity 
analysis  of  built-up  structures.  To  avoid  the  inaccuracies 
of  finite  element  results  on  the  boundary,  the  sensitivity 
information  was  expressed  as  domain  integrals  rather  than 
boundary  integrals.  In  such  an  approach,  design  sensitivity 
formulas  for  each  standard  component  type  (rod,  beam,  plate, 
etc.)  are  obtained,  and  one  can  define  a library  of 
structural  components  that  may  be  assembled  to  form  a built- 
up  structure. 

Another  important  problem  in  shape  optimization  is  the 
determination  of  the  optimal  cross  section  of  a torsional 
shaft.  Banichuk  and  Karihaloo  [22]  present  a solution  for 
an  optimal  cross  section  of  an  isotropic-solid  bar  for 
minimum  weight  and  subject  to  constraints  on  the  torsional 
rigidity  and  bending  stiffness.  Dems  [23]  used  a domain 
perturbation  technique  to  derive  the  optimality  criterion 
for  an  isotropic,  hollow  bar  which  was  designed  for  maximum 
torsional  rigidity  with  a constant  area  and  a fixed  inner 
boundary.  Adali  [24,25]  concluded  that  the  cross  section  of 
an  anisotropic  hollow  bar  with  a high  torsional  rigidity  was 
given  by  a domain  bounded  by  two  geometrically  similar 
ellipses . 

The  torsional  shape  optimization  problem  was  also 
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solved  by  using  an  approach  in  which  the  material 
derivatives  were  employed  to  compute  the  design  sensitivity 
[26] . A finite  element  formulation  was  used  for  the 
analysis.  An  extension  of  this  investigation  [27]  involved 
consideration  of  internal-node  movement  associated  with  the 
regridding  process  on  the  shape  design  sensitivity. 

2.1.2  Boundary  Element  Method  for  Optimum  Shape  Design 

The  first  numerical  treatment  using  the  BEM  in 
structural  shape  optimization  was  proposed  by  Mota  Soares, 
Rodroigues,  Oliveira  Faria  and  Haug  [28].  The  optimization 
of  the  geometry  of  solid  and  hollow  shafts  for  torsional 
requirements  was  formulated  in  terms  of  the  BEM.  Various 
analysis  models,  based  on  constant,  linear  and  quadratic 
boundary  elements  were  used  requiring  only  a discretization 
along  the  boundary.  The  design  objective  was  to  obtain  the 
optimal  cross  section  of  a shaft  with  a given  area,  so  as 
to  maximize  its  torsional  stiffness.  The  discrete  boundary 
element  models  are  much  more  efficient  and  robust  than  the 
corresponding  finite  element  discretizations.  This  is 
largely  attributable  to  the  higher  accuracy  of  analysis 
response,  and  hence  the  more  accurate  sensitivity 
information  obtained.  Furthermore,  there  is  no  need  to 
compute  the  state  variables  in  the  domain,  because  the 
maximum  stress  is  at  the  boundary. 
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Mota  Soares,  Rodrogues  and  Choi  [29,30]  have  also 
developed  the  shape  optimal  design  of  two-dimensional 
elastic  structures  based  on  linear  and  quadratic  boundary 
elements.  The  corresponding  nonlinear  programming  problem 
is  solved  by  a linearization  method.  The  design  objective 
in  this  problem  requires  a minimization  of  the  compliance  of 
the  structure,  subject  to  an  area  constraint.  All  degrees  of 
freedom  of  the  model  are  defined  only  along  the  boundary. 
Thus,  there  is  no  need  for  calculating  displacements  and 
stresses  in  the  domain. 

The  determination  of  the  optimal  shape  for  two-  and 
three-dimensional  elasticity  problems,  based  on  the  boundary 
element  method,  has  recently  been  developed  by  Burczynski 
and  Adamczyk  [31].  A maximal  stiffness  criterion,  together 
with  the  limitation  of  the  volume  of  a body,  have  been  used. 
The  optimality  conditions  were  derived  for  an  optimal 
boundary.  An  iterative  process,  based  on  finite  differences 
and  the  Newton-Raphson  method,  is  used  to  solve  a set  of 
nonlinear  algebraic  equations. 

Choi  and  Kwak  have  developed  a general  method  for  the 
shape  design  sensitivity  analysis  as  applicable  to 
elasticity  problems  using  the  boundary  integral  equation 
formulation  [32].  The  material  derivative  concept  and 
adjoint  variable  method  were  applied  to  obtain  sensitivity 
formulae  for  the  stress  constraint  functions.  Accuracy  of 
the  design  sensitivity  analysis  obtained  from  the  material 
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derivative  approach  was  compared  with  the  conventional 
finite  difference  method.  The  design  of  a fillet  and  an 
elastic  ring  were  used  to  illustrate  the  effectiveness  of 
the  material  derivative  concept.  A smooth  characteristic 
function  was  found  to  be  better  than  a plateau  function  for 
the  localization  of  the  stress  constraint. 

Sandgren  and  Wu  [33]  used  the  BEM  with  substructuring 
concepts  for  two-  and  three-dimensional  structural  shape 
optimization.  Substructuring  in  the  boundary  element  method 
allows  the  shape  changes  of  each  substructure  to  be  isolated 
for  separate  analysis.  The  control  points  on  B-spline 
curves  and  surfaces,  which  were  introduced  to  describe  the 
shape  of  the  domain  were  selected  as  design  variables.  The 
generalized  reduced  gradient  method  was  used  to  solve  the 
nonlinear  programming  problem. 

2 . 2 Formulation  of  Shape  Optimization  Problem 

A general  formulation  of  the  shape  optimization  problem 
usually  involves  minimizing  either  the  structural  mass  or 
stress  by  selecting  shape  design  variables  and 
simultaneously  satisfying  prescribed  geometry  and  response 
constraints  (or  behavior  constraints) . The  mathematical 
formulation  of  a structural  shape  optimization  problem  can 
be  stated  as  follows: 
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Find  xi  i=l , 2 , . . . ,n 

Minimize  F(xjJ 

Subject  to  gj (Xi)  < 0 j =1 , 2 , ,m 

hk(xi)  = 0 k=l,2, . . . ,q 

Xi1  < xi  < xiu  i=l , 2 , . . • , n 

(2.1) 

where  xi  is  a vector  of  shape  design  variables;  F is  the 
objective  function;  gj  and  h^  are  inequality  and  equality 
constraints  on  the  response  quantities;  x^  and  xiu  are  the 
lower  bound  and  upper  bound  of  each  design  variable. 

The  structural  geometry  is  mainly  described  by  the 
definition  of  shape  design  variables.  The  dimensionality  of 
the  design  space  is  dictated  by  the  number  of  design 
variables.  From  the  standpoint  of  the  efficiency  of  the 
optimization  process  and  from  manufacturing  considerations, 
design  variables  are  usually  chosen  as  control  parameters, 
most  typically,  the  coefficients  of  polynomials  or  spline 
curves  in  the  shape  representation. 

The  side  constraints  on  each  design  variable  serve  as 
the  geometry  constraints  which  represent  technological 
limitations  on  the  design  variables  and  reflect  fabrication 
or  analysis  validity  considerations.  The  equality 
constraints  govern  the  structural  response  associated  with 
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the  loading  conditions,  which  are  usually  obtained  from  the 
analysis  tool.  The  inequality  constraints  represent  the 
restrictions  on  structural  response  (eg,  stresses, 
displacements,  or  natural  vibration  frequencies) . In  the 
shape  optimization  problem,  the  von  Mises  stress  function  is 
usually  used  to  calculate  the  equivalent  stress  so  that  each 
boundary  point  has  one  stress  constraint. 

Once  the  optimal  shape  design  problem  is  formulated, 
there  are  two  tasks  that  must  be  considered.  The  first 
relates  to  choosing  the  optimizer  to  achieve  an  optimal 
objective  function  and  satisfy  a set  of  constraints.  The 
second  is  in  the  selection  of  the  tool  to  implement  the 
structural  analysis.  Between  the  optimizer  and  analyzer,  a 
grid  generator  is  used  to  facilitate  data  preparation  for 
the  structure  analysis.  The  grid  generator,  a preprocessor 
of  the  analytical  tool,  converts  the  design  variables 
obtained  from  the  optimizer  to  structural  geometry  variables 
and  generates  the  nodes  for  the  analyzer.  During  the 
optimization  process,  the  optimizer  proposes  trial  design 
variables  to  the  analyzer  through  the  grid  generator,  and 
the  analyzer  calculates  the  response  constraints  and  feeds 
all  analysis  results  back  to  the  optimizer.  Figure  2.1 
illustrates  the  relation  of  the  optimizer,  the  analyzer  and 
grid  generator. 
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Optimizer 


Figure  2 . 1 


The  relation  among  optimizer,  grid  generator, 
and  analyzer. 
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Sensitivity  analysis  deals  with  the  calculation  of  the 
derivatives  of  response  with  respect  to  the  design 
variables.  During  a typical  optimization  seguence,  the 
sensitivity  analysis  accounts  for  most  of  the  CPU  time.  It 
is  precisely  for  this  reason  that  the  issue  of  design 
sensitivity  analysis  in  optimal  shape  design  has  been 
studied  so  extensively.  The  choice  of  proper  sensitivity 
analysis  method  depends  on  the  characteristics  of  both  the 
design  model  and  the  analytical  tool.  A detailed 
introduction  of  sensitivity  analysis  will  be  discussed  in 
Chapter  5 . 


2 . 3 Overview  of  Mathematical  Programming  Method 

The  purpose  of  mathematical  programming  techniques  is  to 
provide  a computer  tool  to  aid  the  designer  in  performing 
various  optimal  design  tasks.  The  use  of  these  techniques 
in  engineering  design  problems  offers  a logical  approach  to 
design  automation.  The  category  of  the  mathematical 
programming  techniques  of  optimization  can  be  classified 
into  two  distinct  types:  (A)  linear  mathematical 

programming,  (B)  nonlinear  mathematical  programming.  In  the 
linear  programming  (LP)  problem,  all  constraints  and  the 
objective  function  are  linear  functions  of  the  design 
variables.  The  simplex  method  is  available  for  solving  the 
linear  programming  problem,  as  is  the  recently  introduced 
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and  still  controversial  approach  proposed  by  Karmarkar  [34]. 
In  these  problems,  the  optimum  always  occurs  on  the  boundary 
of  the  feasible  domain  and  the  exact  global  optimum  is 
reached  in  a finite  number  of  steps.  Although  most 
engineering  problems  are  not  of  the  LP  form,  some  nonlinear 
programming  (NLP)  problems  can  be  approximated  by  linear 
expressions  as  in  the  linear  sequential  programming  (LSP) 
method.  Also  the  search  direction  in  the  NLP  method  often 
requires  the  solution  of  an  auxiliary  LP  problem,  as  in  the 
usable-feasible  search  direction  method. 

In  most  of  the  structural  optimization  problems,  both 
the  objective  function  and  the  constraints  might  be  highly 
nonlinear  functions  or  even  be  implicit  functions  of  the 
design  variables.  The  following  discussions  are  focused  on 
describing  the  characteristics  of  the  nonlinear  optimization 
algorithms  used  in  the  present  work. 

In  the  application  of  NLP  method  in  engineering 
optimization  problems,  convergence  to  the  global  optimum  is 
seldom  assured.  The  design  spaces  in  such  problems 
typically  nonconvex  and  several  local  optima  exist  in  the 
region  of  search.  Another  pitfall  is  that  a poorly  scaled 
design  space  obtained  for  a particular  choice  of  design 
variable  may  result  in  extremely  slow  convergence  of  the 
optimization  algorithm.  One  way  to  increase  the  probability 
of  locating  the  global  optimum  is  to  start  the  optimization 
process  from  several  initial  points. 
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Nonlinear  programming  algorithms  can  be  classified  into 
distinct  categories  of  unconstrained  and  constrained 
optimization  methods.  The  necessary  conditions  for  a local 
minimum  for  unconstrained  problems  is  that  the  gradient  of 
the  objective  function  must  vanish  at  that  point,  and 
furthermore,  the  Hessian  matrix  must  be  positive  definite. 
Here  the  Hessian  matrix  is  defined  as  the  matrix  of  second 
partial  derivatives  of  the  objective  function  with  respect 
to  the  design  variables.  In  a constrained  optimization 
problem,  the  necessary  conditions  for  a design  point  (x*)  to 
be  considered  a local  optimum  are  defined  in  terms  of  the 
Kuhn-Tucker  conditions  of  optimality.  These  are  written  as 
follows : 


x*  is  feasible 
^ j 9j  (x  ) = 0 / A j > 0 

m 1 

VF (x*)  + S Aj  vg-j  (X*)  + s A Vhk(x*)  = 0 
j=l  k=l  K m 

Aj  > 0 , Am+k  unrestricted  in  sign 

where  the  symbols,  j,  k,  m,  gj  and  hk,  are  defined  in 
equations  (2.1);  the  A' s are  the  Lagrange  multipliers  which 
are  uniquely  determined  from  the  gradients  of  the  objective 
function  and  active  constraints;  VF,  Vgj , and  Vhk  are  the 
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first  derivatives  of  objective  function,  inequality 
constraints,  and  equality  constraints,  respectively,  with 
respect  to  the  design  variables. 

The  available  algorithms  for  unconstrained  NLP  can  be 
categorized  on  the  basis  of  the  function  and  function 
derivative  information  required  in  computing  the  search 
direction.  Examples  of  this  category  is  as  follows, 

(1)  Zero-order  methods:  Random  search  method  [35]. 

Powell's  method  [36]. 

(2)  First-order  methods:  Steepest  descent  method  [35]. 

Conjugate  direction  method  [37]. 
Variable  metric  methods  [38]. 

(3)  Second-order  methods:  Newton's  method  [39,40]. 

The  above  methods  each  have  their  own  advantages  and 
drawbacks.  Generally,  the  more  sophisticated  methods 
converge  more  rapidly  and  provide  more  accurate  results. 
However,  the  efforts  involved  in  computing  the  searching 
direction,  particularly  those  related  to  obtaining  gradient 
information,  must  be  taken  into  consideration  when 
determining  the  effectiveness  of  a given  method. 

The  available  algorithms  for  constrained  NLP  are  divided 
into  two  groups.  The  first,  indirect  methods,  converts  a 
constrained  problem  to  an  unconstrained  problem  by  appending 
the  violated  constraints  to  the  objective  by  means  of  a 
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penalty  function.  In  such  an  approach  the  optimum  is 
achieved  by  sequentially  solving  an  unconstrained 
minimization  problem  with  several  values  of  the  penalty 
parameters.  This  approach  is  generally  separated  into  the 
exterior  penalty  function  method,  the  interior  penalty 
function  method,  and  the  extended  penalty  function  method. 
These  methods  are  discussed  in  more  detail  by  Fiacco  and 
McCormick  [41].  The  augmented  Lagrange  multiplier  method 
[42]  utilizes  the  Kuhn-Tucker  conditions  and  includes  the 
Lagrange  multipliers  into  the  exterior  penalty  function, 
improving  the  convergence  of  the  optimization  process. 

The  second  approach  is  a class  of  direct  methods  which 
solve  the  constrained  problem  directly.  The  methods  of 
sequential  linear  programming  (SLP)  method  [43],  feasible 
direction  method  [44],  generalized  reduced  gradient  (GRG) 
method  [45],  and  sequential  quadratic  programming  method 
[46]  belong  to  this  category.  When  the  evaluations  of  the 
objective  function  and  sensitivity  are  costly,  these 
approaches  are  preferable.  The  following  discussion  focuses 
on  those  methods  which  are  used  in  the  present  research. 

The  sequential  linear  programming  method  linearizes  the 
nonlinear  problem  with  a first-order  Taylor  expansion  as 


follows . 
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Minimize  F(x)  = F(x°)  + VF(x°)  • fix 

S.T.  gj (x)  = gj (x°)  + vgj (x°)  • fix  < 0 

hk(x)  = hk(x°)  + Vhk(x°)  • fix  = 0 

Xi1  < Xi+6Xi  < XiU 

where  fix  = x-x°  (2.2) 

The  objective  function  and  the  constraints  are 
linearized  at  the  initial  design  x°.  By  constructing  move 
limits  for  each  design  variable,  the  optimum  of  the 
nonlinear  problem  is  obtained  within  a few  linear 
programming  cycles.  The  method  converges  well  for  fully 
constrained  problems.  However,  if  the  optimum  is  located  at 
a point  where  only  a few  constraints  are  active,  the 
convergence  pattern  sometimes  exhibits  serious  oscillations. 

The  usable  feasible  search  direction  method  belongs  to 
the  category  of  direct  methods.  A FORTRAN  program,  CONMIN 
[47],  based  on  this  method  is  widely  applied  in  engineering 
problems.  In  this  approach,  a direction  vector  {S}  is 
determined  which  satisfies  both  feasible  and  usable 
conditions,  expressed  mathematically  as  {S)T{Vgj}  < 0 and 
{S}t{VF}  < 0,  respectively.  A physical  interpretation  of 
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such  a search  direction  is  that  it  must  point  in  a direction 
which  allows  no  constraint  violation  and  also  result  in  a 
decrease  of  the  objective  function  (for  nonlinear  problems) . 
The  search  direction  vector  { S } is  obtained  as  a solution  of 
the  linear  programming  problem  stated  as  follows: 

Max.  p 

Subject  to  {S}T (Vgj)  + p <0 

{S}T  {VF}  + p < 0 

-{1}  < { S > < (1)  (2.3) 

Here  p is  a positive  scalar  quantity,  larger  value  of  which 
makes  (S)T{VF)  more  negative,  and  forces  the  search 
direction  (S)  to  point  in  a direction  in  which  the  objection 
function  decreases  more  rapidly;  6 j is  a positive  push-off 
factor  which  drives  the  search  direction  away  from  the 
constrained  boundary  and  into  the  feasible  domain. 

The  advantages  of  this  method  are  that  only  a linear 
programming  problem  is  solved  in  determining  the  search 
direction,  and  furthermore,  only  the  active  constraints  need 
to  be  considered  in  calculating  the  gradient  information, 
usually  a time-consuming  process.  However,  this  method  is 
only  well  suited  for  fundamentally  inequality  constrained 
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problems.  If  equality  constraints  are  introduced  into  the 
problems,  this  method  will  not  be  applicable  in  the  form 
presented  above. 

The  generalized  reduced  gradient  method  is  developed  to 
solve  nonlinear  problems  in  the  sense  that  the  simplex 
algorithm  exists  to  solve  the  linear  problem.  This  is 
accomplished  by  adding  slack  variables  to  inequality 
constraints.  The  general  form  of  eqn.  (2.1)  is  thus  stated 
as , 

z n-q  independent  variables 

Find  x = { } 

y m+q  dependent  variables 

Minimize  F(x)  = F(z,y) 

Subject  to  V>k(x)  = 0 k=l , 2 , . . . , m+q 

xi1  < xi  < xiu  i=l , 2 , . . . , n 

xj+n  * 0 j=l,2, . . . ,m  (2.4) 

Differentiating  the  objective  and  constraint  functions  in 
eqn.  (2.4),  the  following  forms  are  obtained. 

dF(x)  = VZF (x)  • dz  + VyF (x)  • dy  (2.5) 

d^k(x)  = VzV>jc(x)-dz  + VyV>k(x)dy  = 0 k=l,m+q 

(2.6) 
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Here,  the  subscripts  z and  y identify  the  gradients  with 
respect  to  the  independent  and  dependent  variables, 
respectively. 

By  solving  for  dy  in  eqn.  (2.6)  and  substituting  dy 
into  eqn.  (2.5),  the  generalized  reduced  gradient,  Vr,  is 
obtained  as  follows. 

dF  (x) 

vrFT(x)  = = VzFt(x)  - vyFT(x)  • [VyV>f  1 • VZV> 

dz 

(2.7) 

The  reduced  gradient  defines  the  rate  of  change  of  the 
objective  function  with  respect  to  the  decision  variable 
with  the  state  variable  adjusted  to  maintain  feasibility. 

The  problem  of  eqn.  (2.4)  is  then  converted  to, 

Min.  f(z)  ; z = [ z± , z2 , . . . , zQ] T 

S.T.  z(1)  < z < z(u)  (2.8) 


where  the  state  variables  y have  been  eliminated  from  the 
original  problem  by  using  the  constraints  ^m(z,y)=0  to  solve 
for  y in  terms  of  z. 

A necessary  condition  for  the  existence  of  a minimum  of 
an  unconstrained  nonlinear  function  is  that  the  elements  of 
the  gradient  vanish. 
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Vrf(x)i  = 


> 0 if  zi  =zj/1^ 
< 0 if  zi  =zj/u^ 
= 0 otherwise 


i 1,2,. ..,Q 

(2.9) 


The  use  of  GRG  method  effectively  solves  the  problems 
with  a few  inequality  constraints  and  several  equality 
constraints.  If  there  are  a large  number  of  equality 
constraints,  the  computer  storage  becomes  prohibitive. 


CHAPTER  3 


BOUNDARY  ELEMENT  METHODS  FOR  STRUCTURAL  ANALYSIS 

3 . 1 Introduction 

The  boundary  element  method  is  an  alternative  approach 
to  obtaining  a solution  to  a given  set  of  governing 
equations  in  continuum  mechanics.  The  technique  consists  of 
the  transformation  of  the  partial  differential  equation 
describing  the  behavior  of  the  unknowns  inside  and  on  the 
boundary  of  the  domain  into  an  integral  equation  relating 
only  boundary  values.  The  geometry  of  the  boundary  is  then 
discretized  into  several  elements  to  formulate  a set  of 
linear  system  equations.  Numerical  solutions  of  these 
equations  with  prescribed  boundary  values  allow  a 
computation  of  the  unknown  responses  along  the  boundary.  If 
values  at  the  internal  points  are  required,  they  are 
calculated  from  the  boundary  information. 

Since  all  numerical  approximations  take  place  only  at 
the  boundaries,  the  dimensionality  of  the  resulting  linear 
system  of  equations  is  considerably  smaller  than  that 
obtained  in  the  finite  element  method.  Hence,  the  amount  of 
computation  required  is  substantially  reduced,  as  is  the 
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amount  of  labor  involved  in  defining  the  geometry  of  the 
problem.  These  favorable  characteristics  contribute  towards 
the  recognition  of  the  boundary  element  method  as  an 
efficient  analysis  tool  for  the  optimization  process. 

The  interest  in  engineering  applications  of  the  BEM 
increased  significantly  in  the  latter  part  of  the  last 
decade,  and  there  have  been  several  recent  contributions  to 
the  literature  in  this  field  [48].  In  this  chapter,  a brief 
introduction  and  application  of  the  boundary  element  method 
is  introduced.  Two-dimensional  elastic,  torsional  shaft 
problems,  including  singly-  and  doubly-connected  prismatic 
bars,  and  a small-strain  elastostatics  problem  is  solved  by 
the  BEM,  and  these  implementations  are  presented  here  for 
completeness . 


The  BEM  for  the  potential  problem  is  based  on  Green's 
formula,  eqn.(3.1),  which  allows  the  formulation  of  certain 
boundary  value  problems  as  integral  equations  [48]. 


3 • 2 Application  in  Torsional  Shaft  Problems 


(w- 


du 


• u)  dr 


(3.1) 


Since  the  governing  equation  of  the  torsional  problem  is  an 
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ellip"tical  partial  differential  equation,  it  may  be  treated 
as  a potential  problem. 

Let  u denote  the  stress  function  of  a torsional 
elastic  problem.  This  stress  function  should  satisfy  the 
Poisson's  equation  as  follows: 

Lu  = vzu  = -b  in  0 (3.2) 

Here  L is  the  differential  operator,  V2  is  the  Laplace 
operator,  and  b is  a prescribed  constant  ( b=2  for  the 
problem  of  torsion  of  an  elastic  bar) . 

The  boundary  conditions  corresponding  to  this  torsional 
problem  are  of  two  types,  as  shown  in  Figure  3.1.  These 
types  are 

(a)  essential  condition,  u = u on  r^,  where  u is  a 
prescribed  value  of  stress  function; 

(b)  natural  condition,  q = — = q on  T2,  and  q is  a 

prescribed  value  of  normal  derivative  over  the  boundary. 

The  total  boundary  is  denoted  by  = r 2 + T3« 

Assume  that  w is  a fundamental  solution  of 

Lw  = -6  (x,xQ)  (3.3) 

where  S(x,x0)  is  the  Dirac  delta  function  which  has  the 
property 
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Figure  3.1  Definition  of  the  boundary  element  method  for 
potential  type  problem. 
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Jn  u(x)  8 (x,xQ)  dx 


u(xQ)  , xD  e o 
1/2  u(xQ)  , xQ  e r 


(3.4) 


The  fundamental  solution  of  the  above  problem  is  of  the 
form, 


w = i 


1 


4 r 
-1 

2n 


for  3-D 
for  2-D 


(3.5) 


where  r=  | x - xD | . 

Applying  Green's  identity,  eqn.  (3.1),  and  making  use 
of  the  fundamental  solution  in  eqn.  (3.4)  allows  one  to 
write  an  expression  for  u at  a point  Xp  as  follows. 

f 3u  3w  f 

u(xp)  = Jr  (w--^ — -u)  dr  + b-w  dQ  (3.6) 


It  is  observed  that  the  fundamental  solution,  eqn. 
(3.5),  is  valid  for  any  point  in  the  domain  n.  However,  if 
xQ  approaches  the  boundary  r , then  the  solution  at  xQ 
becomes  singular # In  order  to  formulate  the  problem  by  the 
BEM  approach,  the  singularity  must  be  specially  handled  by 
mathematical  analysis.  For  instance,  in  a 3-D  problem,  if 
the  boundary  is  assumed  to  be  smooth,  one  may  cut  out  a 
small  hemispherical  region  of  the  boundary  to  analyze  the 
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singular  point  (see  Figure  3.2).  The  boundary  integration 
around  the  singular  point  Xp  is  calculated  as 


r aw  f i 

iim  ( Je  U-—  dr  ) = lira  ( -L  u — — - dr  ) = -1/2  u 

e >0  °n  e >0  4,re 


lira  ( L 
€ >0  J 


i au 

dr  ) 

4?rr  an 


0 


where 

dr  = r sin  6 ■ dc/>  ■ d?r 


Thus,  if  Xp  is  on  the  boundary,  equation  (3.6)  becomes 


f aw  r au  r 

1/2 • u (xp)  + Jr  u dr  = r -w  dr  + r b-w  do 

an  an 

(3.7) 


To  implement  the  BEM  on  a digital  computer,  the  boundary 
of  the  structure  geometry  is  discretized  into  n segments. 
Equation  (3.7)  can  then  be  written  in  a discrete  form  as, 


n n 

l/2-Ui(Xp)  + 2 Uj  • Hij  = 2 qj  • Gij  + Bi  (3.8) 

j=l  j=l 


where, 

f aw 

Hij  = Jr  dr 

an 


G^j  — J*  r w dr 
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Figure  3.2  Boundary  surface  assumed  to  be  hemispherical  for 
integration  purpose. 
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Bi  = Jn  b-w  do 

The  first  two  terms  of  the  left  hand  side  in  equation 
(3.8)  are  combined  in  one  term,  still  represented  by  H^j  . 
The  system  equation  can  now  be  written  as, 

n n 

S Hij  • Uj  = S Gij  • qj  + Bi  (3.9) 

j=l  j=l 

or  can  be  expressed  in  matrix  form  as  follows. 

[H] • (U)  = [G] • (Q)  + { B)  (3.10) 

The  boundary  conditions  described  above  are  only 
applicable  to  the  singly-connected  torsional  bar.  When  the 
problem  involves  torsion  of  a prismatic  bar  with  a doubly- 
connected  domain  (Figure  3.3),  the  boundary  conditions  will 
change  as  follows  [23], 

u = 0 on  r 

u = uQ  on  rQ 

3u 

ds  = 2 Do  on  rQ  (3.11) 

3 n 

where  D0  is  the  area  enclosed  within  the  interior  boundary, 
and  K,  the  torsional  rigidity,  is  expressed  as 
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Figure  3 . 3 


A multiply-connected  domain. 
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K = 2 


I 


n u dft  + 2 uQ  • D0 


(3.12) 


For  the  multiply-connected  domain,  one  unknown  constant 
uD  and  one  additional  equation  (3.11)  are  introduced  into 
the  linear  equation  (3.10).  The  dimensions  of  the  matrices, 
H and  Q,  are  then  modified  as  (n+l)*(n+l),  and  the 
dimensions  of  the  matrices  U,  Q,  and  B are  changed  to 
(n+1) *1. 


The  boundary  element  method  for  elasticity  problems  is 
based  on  Somigliana's  identity  [3],  and  is  given  by. 


where  uj (£ ) is  the  displacement  at  point  £ in  j direction. 
The  quantities  pj (x)  and  bj (x)  are  the  actual  state  of  the 
traction  and  the  body  force  at  point  x,  respectively; 
Gij(£/X)  and  F^(£,x)  represent  the  Kelvin  solution  for 
displacement  and  traction  at  £ in  the  j direction  due  to  a 


3 . 3 Application  in  Elasticity  Problems 


(3.13) 
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unit  force  applied  at  x in  the  i direction.  These  functions 
can  be  expressed  as  follows, 


1 

8ttG(1-j/) 


t (3-4./) 


ln(— £)  -fiij  + 


dr  dr 

ax.  dXjJ 


(3.14) 


Fij<*'x>  = 


4?r  ( 1 -v  ) r L an 


r ar  ,,,,,,  , . ar  dr 

t — ((i-2-)' ‘ij  + 2 ) 


,,  _ . ,3r  „ ar  , , 

(1-2l/) ' (ax.  'nj  ax.  'ni}  ] 

1 D 


(3.15) 


where, 

r = | x - e I - 

If  the  field  point  £ is  close  to  the  boundary,  one  arrives 
at  an  expression  involving  only  the  boundary  unknowns. 
Since  the  singularity  point  is  introduced  during  the 
boundary  integration,  the  integral  equation  (3.13)  assumes 
the  form, 


uj  U) 


Jr  [ P(x)  Gij (e,x)  - Fij (£,x)  uj (x) ] dr (x) 


b-Gij  (£,x) 


do  (x) 


(3.16) 


where  C. . 

ID 


can  be  expressed  as  follows, 
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C.  . = S . . - I . . 
13  ID  ID 


(3.17) 


Here,  1^^  is  defined  in  terms  of  9 ]_  and  9 2 (Figure  3.4),  and 


is  written  as, 


1 

8tt  ( l~i/  ) 


4 (1-u)  (n  + 9 2 -9  i ) 

+ sin20x  - sin20z 
COS20Z  - 0032^ 


cos20z  “ COS20-L 

4 (1-u) (n+9 2 ~9 1 ) 

+sin20z  “ sin20x 


where  C-^j  is  the  coefficient  that  depends  on  the  geometry  of 
the  boundary  at  point  £,  i.e.,  equation  (3.16),  or  may  be 
determined  from  rigid  body  motion  [48]. 

The  numerical  solution  of  eqn.  (3.13)  is  facilitated  by 
discretizing  the  boundary  into  N elements.  The  values  of 
displacement  u and  traction  p over  each  element  are 
approximated  by  using  the  interpolation  function  <t> , 
expressed  in  terms  of  the  nodal  values  un  and  pn  as  follows. 

u = <t>T  un  (3.18) 


P = 


n 


P 


(3.19) 


Substitution  of  expressions  for  u and  p into  equations 
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Figure  3.4 


Definition  of  domain  boundary  orientation  as 
defined  by  angles  8±  and  82- 
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(3.16)  and  summing  over  all  the  N segments,  results  in  an 
alternative  matrix  notation  as  follows. 


where  [H]  and  [G]  contain  the  results  of  all  the 
integrations.  These  integrals  may  be  evaluated  by  any 
suitable  method  such  as  the  Gauss  quadrature  integration 
scheme.  The  system  of  equations  (3.20)  can  be  rearranged  in 
such  a way  that  all  unknowns  are  written  on  the  left  hand 
side  in  a vector  of  unknowns  (x),  such  that 


where  [A]  is  an  N*N  matrix,  and  (C)  is  a vector  of 
prescribed  boundary  values. 

Once  the  boundary  values  are  calculated,  one  can  obtain 
stresses  at  any  internal  point  by  using  the  stress  and 
displacement  relation. 


[H]  • (U>  = [G]  • (P)  + (B) 


(3.20) 


[A]  (x)  = (C) 


(3.21) 


(3.22) 


where  A and  G are  Lame'  constants.  Differentiating  equation 
(3.13)  and  substitution  into  the  equation  (3.22),  the 
stresses  at  the  internal  points  are  obtained  as  follows. 
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2Gi/ 

"ij  = Jr  ‘ TX 


^lk  3Gik  3Gik 

6..  — + G ( — + — 1 — )}  Pk  dr 

J a xi  sx-i  axi 


r 2Gj/ 

{ 

JQ  1-2is 


aGlk  aGik  aGik 

S-M  — ~ + G ( — + — 1 — ) } bk  dQ 

J a xi  axj  axi 


I 


2Gv  aF  3F.  aF. 

{ 6.,  + G( — + — 1 — ) } Uk  dr 

c i-2i^  J a xi  axj  axj 


(3.23) 


In  general,  the  stress  tensor  on  the  boundary  is  also 
important  for  structural  shape  design  problem.  The 
derivation  of  the  stress  tensor  along  the  boundary  is  shown 
in  the  Appendix. 


3 . 4 Numerical  Results 


Numerical  implementations  of  the  boundary  element  method 
are  illustrated  by  an  elliptic  torsional  shaft  and  a 
cylindrical  shell  in  a state  of  plane  strain.  Since  the 
performance  of  structural  analysis  may  affect  convergence 
during  the  shape  optimization,  both  numerical  results  are 
compared  with  analytical  solutions.  Computer  programs  for 
the  solution  of  the  isotropic  torsional  problem  and  two- 
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dimensional  elastostatics  problem  are  generated  in  this 
study.  These  programs  are  available  for  both  constant  and 
linear  elements  along  the  boundary. 

As  mentioned  earlier,  the  program  requires  that  only 
the  boundary  of  the  domain  be  discretized.  There  is  no 
requirement  of  a special  assembler  subroutine  as  in  the 
finite  element  method.  Furthermore,  not  only  is  the  data 
preparation  much  simpler  but  also  the  order  of  the  resulting 
linear  system  equations  is  considerably  lower  than  that  in 
the  finite  element  approach.  However,  the  sparseness  and 
the  symmetry  of  the  influence  matrix  [A]  are  generally  lost. 

Figure  3.5  (a)  shows  the  solution  of  the  stress 
function  for  a torsional  prismatic  shaft  of  an  elliptic 
cross  section.  Figure  3.5  (b)  displays  the  stress 
distribution  along  the  elliptic  boundary.  The  comparison  of 
boundary  stresses  between  computational  and  analytical 
solutions  are  listed  in  table  3.1. 


Table  3 . 1 Comparison  of  stresses  along  the  boundary 


X 

y 

r/G/3 

(analytical) 

r/G/3 

( computational ) 

error 

% 

5 

0 

2.6471 

2.6493 

0.083 

4 . 05 

1.76 

3 .3632 

3 . 3749 

0.347 

2 . 5 

2 . 6 

4 . 0435 

4.0583 

0.366 

0.52 

2 . 98 

4 .3964 

4 .4116 

0.346 
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Figure  3.5  (a)  Solution  of  stress  function  for  elliptic 
torsional  shaft. 


stress 


Figure  3.5  (b)  Shear  stress  distribution  along  the 
boundary. 
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The  numerical  results  of  a prismatic  bar  under  torsion 
with  a doubly-connected  domain,  is  also  examined  in  this 
section.  Assume  that  A and  B define  the  major  and  minor 
lengths  of  the  outer  elliptic  boundary,  respectively,  and 
the  inner  elliptic  cutout  is  represented  by  the  major  axis  a 
and  a minor  axis  b.  Two  examples,  figure  3.6,  for  different 
major  and  minor  lengths  are  used  to  illustrate  the  numerical 
results.  In  the  tables,  uQ  represents  the  inner  unknown 
constant  in  equation  (3.11),  K is  torsional  rigidity  as 
expressed  in  equation  (3.12),  and  the  maximum  and  minimum 
stresses  along  the  boundary  are  also  listed.  The  torsional 
rigidities  for  this  case  are  used  to  compare  the  numerical 
solutions  with  analytical  solutions. 

Figure  3.7  shows  a cylindrical  shell  subjected  to  the 
internal  pressure  of  Pa  = 100  units  and  external  pressure  of 
Pb  = “25  units  acting  on  the  surfaces  at  ra  = 3 and  rb  = 8 , 
respectively.  The  problem  can  be  assumed  to  be  one  of  plane 
strain.  Variations  in  both  the  computational  and  analytical 
solutions  of  and  t along  the  radial  direction  are 

shown  in  figure  3.8  (a)  and  (b) , respectively.  It  is 
apparent  that  the  accuracy  of  the  numerically  computed 
stress  at  internal  points  near  the  boundary  is  lost.  Since 
the  internal  stress  is  determined  from  the  derivatives  of 
the  integral  representation  for  displacement,  (equation 
(3.23)),  the  term  containing  r’ 2 obtained  from  the 
derivation  of  Equation  (3.23)  in  the  integral  stress 


stress 
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Figure  3.6  Stress  distributation  along  the  multiply- 

connected  elliptic  domain  of  the  torsioal  shaft. 
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Figure  3 . 6 Continued 
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Figure  3. 


P 


b 


Description  of  cylindrical  shell  subjected  to 
internal  pressure  Pa  and  external  pressure  Pfc,. 
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Figure  3.8  (a)  Stress  distribution  r „ versus 

cylindrical  shell.  xx 
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Figure  3.8  (b)  Stress  distribution  r versus 

cylindrical  shell.  ^ 
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results  in  a singularity  for  small  values  of  r.  Due  to  such 
a strong  singularity,  integration  points  cannot  be  allowed 
to  approach  the  boundary.  Stress  at  such  a point  can  be 
obtained  from  the  stress  computed  either  by  using  the 
interpolation  method  or  by  decreasing  the  singularity  by  one 
order  [ 49 ] . 


CHAPTER  4 


BOUNDARY  ELEMENT  METHOD  FOR  OPTIMUM  SHAPE  SYNTHESIS 


4 . 1 Selection  of  Shape  Design  Variables 

The  choice  of  shape  design  variables  has  a rather 
critical  influence  in  determining  the  character  of  the 
boundary  in  a shape  optimization  problem.  A proper 
selection  of  shape  design  variables  not  only  substantially 
influences  the  characteristics  of  the  design  space  in  the 
nonlinear  programming,  but  also  avoids  undesired  or 
impractical  shape  design.  Figure  4.1  shows  the  final  shape 
obtained  in  a the  fillet  design  problem  involving  the  use 
of  a boundary  element  method.  The  jagged  shape  mainly 
results  from  improper  grid  distribution  along  the  boundary. 
This  often  happens  when  coordinates  of  the  boundary  nodes 
are  defined  as  design  variables.  Thus,  either  from  the 
standpoint  of  manufacturing  and  practicality  of  a design,  or 
from  the  view  of  computational  savings,  the  selection  of 
design  shape  variables  plays  a very  important  role  in  the 
shape  optimization  problem. 

Several  representations  of  shape  design  variables  have 
been  proposed  to  describe  the  structural  boundary  in  a 
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Design  generated  when  node  locations  are  used  as 
design  variables 


Figure  4 . 1 
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design  domain.  Except  for  the  use  of  node  coordinates  to 
represent  the  boundary,  various  shape  representations  of 
structural  boundaries  have  been  summarized  in  reference  [1]. 

The  polynomial  representation  has  been  commonly  used  to 
represent  boundaries  in  optimum  shape  design.  After  taking 
the  coefficients  of  the  polynomial  as  design  variables,  the 
design  shape  is  then  characterized  by  a set  of  polynomials 
[50,51].  Dems  [23]  describes  the  design  shape  by  a set  of 
prescribed  shape  functions  and  a set  of  shape  parameters. 

The  optimization  procedure  is  reduced  to  the  determination 
of  these  parameters.  Since  the  use  of  high-order 
polynomials  usually  results  in  an  oscillatory  shape,  spline 
representations  composed  of  several  low-order  elements  which 
are  combined  to  maximize  smoothness  are  used  to  eliminate 
this  problem  [49].  Results  obtained  in  [52]  also 
demonstrate  the  effectiveness  of  spline  representations  in 
producing  better  accuracy  in  the  sensitivity  information  in 
comparison  to  an  approach  in  which  piecewise  liner 
representation  of  the  boundary  is  used. 

The  design  element  concept  was  first  introduced  by  Imam 
[16].  The  design  elements  are  described  by  a set  of  key 
nodes  that  control  the  geometry.  Braibant  and  Fleury  [53] 
used  Bezier  and  B-spine  blending  functions  to  represent  the 
boundary  of  design  elements. 

Many  selections  of  shape  design  variables  are  described 
above.  However,  the  literature  does  not  elaborate  on  how  to 
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choose  the  order  of  the  polynomial  when  such  a polynomial 
representation  is  used.  Furthermore,  in  some  problems  it 
may  be  advantageous  to  remove  material  from  within  the 
domain.  In  such  situations,  it  is  not  always  obvious  what 
kind  of  shape  is  preferred,  what  shape  such  an  interior 
cutout  should  have,  or  where  such  a cutout  should  be 
located.  To  obtain  a better  understanding  of  these 
problems,  a guideline  for  choosing  proper  shape  design 
variables  based  on  the  domain  stress  contours  is  studied  in 
this  section.  This  study  is  limited  to  a class  of  design 
problems  where  the  analysis  domain  is  linear  elastostatics . 

At  the  very  outset  of  generating  a shape  design  modeling 
process,  the  boundary  and  loading  conditions  are  generally 
specified.  Assuming  an  initial  design  domain  that  satisfies 
both  conditions,  the  von  Mises  equivalent  stress  in  this 
domain  can  be  calculated  by  any  analytical  tool.  In  the 
present  study,  the  boundary  element  method  was  used  for  this 
task. 

An  equivalent  stress  contour  plot  is  a useful 
representation  to  gain  a better  understanding  of  the  load 
distribution  in  the  structural  domain.  By  connecting  the 
equal  equivalent  stress  ratio  points  in  the  structural 
domain,  where  the  stress  ratio  is  the  von  Mises  stress 
divided  by  the  allowable  stress,  the  stress  distribution  on 
the  structural  domain  can  be  easily  gauged,  as  shown  in 
figure  4.2.  If  each  equivalent  stress  ratio  line  is 
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distribution  in  the 
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considered  to  be  a stress  path,  it  is  apparent  that  a small 
amount  of  stress  flow  passes  through  the  low  stress  regions. 
When  reviewing  the  problem  from  the  standpoint  of  optimum 
shape  design,  the  low  stress  regions  may  justifiably  be 
treated  as  regions  from  which  material  can  be  removed. 
However,  there  still  exist  problems,  including  those  related 
to  the  existence  of  low  stress  regions,  and  the 
impracticality  of  the  boundary  shapes  if  all  such  regions 
were  removed  in  the  shape  design  problem. 

As  discussed  earlier,  the  eguivalent  stress  contour  plot 
provides  valuable  information  for  optimum  shape  design.  To 
use  the  information  effectively,  several  basic  shape 
elements,  such  as  line  elements,  circle  elements,  triangular 
elements  and  parabolic  elements,  were  defined  and  made 
available  for  optimum  shape  design.  Each  element  was 
handled  by  control  parameters,  such  as  locations  of  the  two 
end  points  for  a line  element,  the  radius  or  the  location  of 
the  center  point  for  a circle,  or  one  end  point  slope  and 
the  location  of  the  two  end  points  for  a parabolic  element. 
Guided  by  the  equivalent  stress  distribution,  these  basic 
elements  were  introduced  into  the  low  stress  regions,  and 
new  structural  boundaries  defined  by  these  elements.  The 
sizes  of  these  elements  were  determined  by  the  control 
parameters.  In  other  words,  the  control  parameters  were 
considered  to  be  design  variables  during  the  shape 
optimization  procedure. 
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An  example  of  stress  contour  based  representations  of 
shape  design  boundaries  will  be  shown  in  Chapter  6,  where 
this  approach  was  combined  with  grid  adaptation  and  grid 
refinement  techniques.  Depending  on  the  structural  geometry 
and  the  loading  conditions,  the  numerical  results  of  the 
boundary  element  method  are  seriously  affected  by  grid 
distribution.  These  problems  are  most  frequently 
encountered  in  structures  where  regions  of  stress 
concentration  exist. 

In  the  following  section,  several  shape  design  examples 
are  presented  using  an  approach  that  combined  nonlinear 
programming  optimization  methods  with  the  boundary  element 
method.  These  shape  design  problems  include  torsional 
shafts  and  structural  components  that  are  characterized  by  a 
state  of  plane  stress,  and  for  which  previous  finite  element 
based  solutions  are  available  in  the  literature.  Design 
variables  chosen  in  this  study  are  nodal  coordinates  and 
coefficients  of  low-order  polynomials.  Since  the  grid 
adaptation  and  grid  refinement  techniques  are  not  applied  in 
these  specific  examples,  the  grids  used  in  two-dimensional 
elasticity  problems  are  evenly  and  densely  located  along  the 
structural  boundaries. 

4 . 2 Optimum  Shape  Design  for  Torsional  Shaft 


The  first  example  in  this  sequence  of  optimum  designs 
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involves  the  determination  of  the  boundary  of  a prismatic 
bar  for  maximum  torsional  rigidity,  subject  to  a constraint 
of  a fixed  cross  sectional  area  of  50  units.  An  additional 
requirement  that  the  moment  of  inertia  about  the  horizontal 
axis  be  greater  than  350  units  was  also  investigated  in  this 
example.  The  definitions  of  design  variables,  geometry 
descriptions,  and  the  initial  design  for  this  shaft  are 
shown  in  figure  4.3.  Six  design  variables  are  used  to 
obtain  the  optimum  design,  using  both  a linear  piecewise 
approximation  to  the  constraints  and  objective  function,  and 
a fully  nonlinear  treatment  of  these  functions.  The  latter 
approach  directly  converges  to  the  optimum  design,  while  the 
former  approach  results  in  oscillation  near  the  optimum 
point.  The  final  optimum  shapes  without  and  with  the 
constraint  on  moment  of  inertia  about  a horizontal  axis,  are 
shown  in  figures  4.4  and  4.5,  respectively.  The  optimum 
values  of  all  design  variables,  the  torsional  rigidity, 
cross  sectional  area,  and  bending  inertia  about  the 
horizontal  axis  for  the  initial  and  final  designs  are  also 
shown  in  each  figure. 

The  second  example  is  an  extension  of  the  previous 
problem,  with  an  interior  hole  introduced  in  the  torsional 
shaft.  This  is  an  example  of  a multiply-connected  domain. 
The  initial  radii  of  outer  and  inner  circle  boundaries  are  8 
and  3 units,  respectively.  The  optimization  problem  is 
simplified  to  define  fourteen  design  variables  in  the  first 
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Design  variables 

X1  x2  x3  x4  x5  x6 

6 6.26  4.97  4.97  6.36  6 

Area  = 85.26 

Bending  inertia  on  horizontal  axis  = 682.30 
Torsional  rigidity  = 567.14 


Figure  4.3  Initial  design  and  design  variable  definition 
for  torsional  shaft 
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Design  variables 

X1  x2  x3  x4  x5  x6 

5.45  4.88  4.71  3.36  3.24  3.26 

Area  = 50.03 

Bending  inertia  on  horizontal  axis  = 350.81 
Torsional  rigidity  = 306.01 


Figure  4.4  Final  design  for  torsional  shaft  with  bending 
inertia  constraint. 
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Design  variables 

X1  x2  x3  x4  x5 

4.43  4.21  4.04  4.04  4.21 

Area  = 50.10 

Bending  inertia  on  horizontal  axis 
Torsional  rigidity  = 325.54 


x6 

4.23 

= 206.98 


Figure  4.5  Final  design  for  torsional  shaft  without  bending 
inertia  constraint. 
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Figure  4.6  Design  variable  definition  for  multiply- 
connected  torsional  shaft. 
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quadrant,  as  shown  in  figure  4.6.  The  torsional  rigidity 
and  bending  inertia  requirements  were  similar  to  the  solid 
shaft  with  additional  requirements  to  constrain  the  total 
and  inner  hole  section  areas  to  40  and  20  units, 
respectively.  A total  of  48  nodes  were  used  in  the 
analysis,  with  converged  results  obtained  in  10  iterations. 
The  initial  and  final  designs  are  shown  in  figure  4.7. 

4 • 3 Optimum  Shape  Design  in  2-D  Elasticity  Problems 

Two  plane  stress  problems  are  attempted  in  this 
section.  A torque  arm  shown  in  figure  4.8,  was  sized  for 
minimum  weight  and  a maximum  allowable  von  Mises  stress  of 
620  Mpa . The  maximum  stresses  in  this  structure  occur  at 
the  boundary.  A total  of  133  linear  boundary  elements  were 
used  in  the  analysis,  with  six  design  variables  used  to 
define  the  geometry.  Three  additional  geometry  constraints 
were  imposed  to  prevent  boundary  intersections  during 
redesign.  The  optimization  problem  was  solved  using  an 
exterior  penalty  function  approach,  with  the  evolution  of 
the  design  shown  in  figure  4.9. 

The  second  example  in  this  class  of  plane  stress 
problems  deals  with  the  minimum  area  sizing  of  a fillet  to 
sustain  a uniaxial  load  shown  in  figures  4.10  and  4.11.  Two 
different  sets  of  design  variables  were  used.  The  first  set 
of  results  (figure  4.10)  were  obtained  for  a single  shape 
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Figure  4.7  Initial  and  final  design  for  multiply-connected 
torsional  shaft. 
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Figure  4.8  Design  variable  definition  and  loading 
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Figure  4.9  Iteration  history  for  torque  arm. 
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Figure  4.10 


Initial  and  final  design  of  fillet  (parabolic 
interpolation) 


Figure  4.11  Initial  and  final  design  of  fillet  (elliptic 
interpolation) 
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variable  which  defines  the  location  at  which  a parabola  from 
point  A would  meet  the  horizontal  line  OB  with  zero  slope. 
The  second  set  of  results  shown  in  figure  4.11  was  obtained 
by  determining  the  location  of  points  A and  B on  the 
vertical  and  horizontal  axis  that  could  be  connected  by  an 
elliptical  curve,  and  satisfy  the  von  Mises  stress 
requirements.  This  problem  is  also  attempted  by  defining 
design  variables  as  radial  locations  from  a fixed  reference 
point  figure  4.1.  An  optimum  design  could  not  be  obtained 
for  this  set  because  of  severe  oscillations  in  the  domain 
boundary.  This  problem  also  reinforces  the  need  to  use  an 
adaptive  node  generation  capability  in  the  BEM  based  optimum 
design  strategy. 


4 . 4 Discussions 


In  this  chapter,  we  have  successfully  introduced  the 
boundary  element  method  in  the  shape  optimization  problem. 

A preliminary  introduction  to  the  definition  of  interior 
cutouts  in  planar  stress  problems,  which  is  based  on  stress 
distribution  in  the  domain,  is  also  presented.  The  use  of 
basic  boundary  shapes  is  advocated  in  such  problems.  The 
applications  of  using  the  boundary  element  method  for  shape 
optimal  design  has  been  demonstrated  in  several  numerical 
examples . 

The  gradient  information  used  in  the  optimization 
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process  for  each  of  the  problems  was  calculated  by  a finite 
difference  approach.  The  sensitivity  analysis  consumes  most 
of  the  CPU  time,  especially  when  the  finite  difference 
method  is  used.  Seeking  a more  efficient  approach  for 
sensitivity  computations  becomes  an  important  task  in  the 
shape  optimization  problem.  An  alternative  approach  of 
utilizing  the  characteristic  of  the  boundary  integral 
equation  to  obtain  the  gradient  information,  is  studied  in 
Chapter  5. 

The  grids  used  in  the  two-dimensional  elasticity 
problems  in  this  chapter,  were  densely  and  evenly  generated 
along  structural  boundaries.  Therefore,  the  order  of  the 
system  of  linear  equations  obtained  from  the  boundary 
element  method  becomes  very  large.  Obtaining  comparable 
levels  of  accuracy  by  reducing  the  number  of  grid  points 
using  the  boundary  element  approach  is  extremely  desirable. 
This  would  yield  substantial  savings  in  the  time  required 
the  solution  of  the  system  equations,  a process  that  is 
repeated  extensively  in  both  analysis  and  sensitivity 
computations.  Toward  this  end,  a grid  refinement  and 
adaptation  technique  is  highly  recommended  when  using  the 
boundary  element  method,  and  will  be  discussed  in  Chapter  6. 


CHAPTER  5 


SENSITIVITY  ANALYSIS  FOR  BOUNDARY  ELEMENT  APPLICATIONS 

5 . 1 An  Overview  of  Sensitivity  Analysis  Methods 

In  most  structure  shape  design  processes,  the  design 
sensitivity  analysis  comprises  an  important  component. 

Design  sensitivity  analysis,  which  involves  the  calculation 
of  the  derivatives  of  the  objective  function  and  constraints 
with  respect  to  the  design  variables,  provides  essential 
information  required  to  couple  mathematical  optimization 
with  structural  analysis  procedures.  Since  the  calculation 
of  sensitivity  information  is  computationally  demanding,  the 
choice  of  a suitable  approach  for  sensitivity  analysis 
becomes  very  important.  There  are  several  approaches 
available  to  perform  shape  design  sensitivity  by  use  of  the 
finite  element  method,  but  there  is  only  limited  experience 
in  the  use  of  the  boundary  element  method  for  this  purpose. 
Sensitivity  analysis  for  the  optimum  shape  design  problem 
involving  the  use  of  the  finite  element  method  for  analysis, 
is  generally  classified  into  four  groups:  (1)  the  direct 

difference  approach,  (2)  the  implicit  differentiation 
approach,  (3)  the  finite  dimensional  approach,  and  (4) 
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the  continuum  derivatives  approach. 

(1)  Direct  difference  approach  method 

The  direct  difference  method  is  commonly  used  in  member 
sizing  optimization.  Botkin  [14]  used  this  approach  to 
approximate  the  shape  design  sensitivity.  By  analyzing  the 
structure  for  both  old  design  variables,  x,  and  new  design 
variables,  x + Ax,  with  given  perturbation  Ax,  the  design 
gradient  is  then  calculated  by  the  finite  difference. 


The  procedure  of  this  method  is  straightforward,  but  is 
costly  when  the  number  of  design  variables  becomes  large, 
particularly  if  the  calculations  of  the  objective  function 
and  the  constraints  are  expensive.  The  accuracy  of  the 
derivatives  is  dependent  upon  the  step  size  of  perturbation. 

(2)  Implicit  differentiation  method 

This  approach  is  usually  applied  in  the  finite  element 
method  [54].  The  finite  element  equilibrium  equation  is 
written  as 


df  (x) 


f ( X+AX)  - f(x) 


(5.1) 


dx 


AX 


[K]  {u}  = {F} 


(5.2) 


where  [K]  is  stiffness  matrix,  (u)  is  displacement  vector, 
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and  {F}  is  the  external  load  vector. 

Differentiating  eqn.  (5.2)  with  respect  to  a design 
variable,  x,  yields  the  following  expression: 

3u  3F  <3K 

[K]  { } = { } - [ ] {u}  (5.3) 

ax  ax  ax 

The  derivative  (au/ax)  can  then  be  used  to  obtain  the 

sensitivity  of  the  stress  response.  The  use  of  such  an 

equation  circumvents  the  necessity  of  solving  for 
displacements  for  perturbed  and  unperturbed  values  of  the 
design  variables.  The  calculation  of  [aK/ax] , however,  is 
expensive,  especially  if  changes  in  the  boundary  may  result 
in  large  changes  in  the  stiffness  matrix.  Furthermore,  the 
changes  in  shape  may  result  in  a distortion  of  the  finite 
element  grid,  and  the  accuracy  of  response  solution  and 
hence  also  the  sensitivity  analysis,  may  deteriorate 
significantly. 

(3)  Finite  dimensional  approach 

If  the  coordinates  of  nodes  are  considered  to  be  the 
design  variables  of  the  shape  optimization  problem,  then  the 
finite  element  stiffness  matrix  is  clearly  dependent  on  the 
design  variables.  One  can  evaluate  the  derivatives  of  the 
stiffness  matrix  and  the  load  vector  with  respect  to  the 
shape  design  variables,  and  the  derivatives  of  the 
constraints  can  then  be  obtained  by  the  adjoint  variable 
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technique  [8,9,55].  The  number  of  adjoint  equations  is  the 
same  as  the  number  of  constraints  and  is  independent  of  the 
number  of  design  variables.  Hence,  if  the  number  of 
constraints  is  less  than  the  number  design  variables,  this 
approach  has  an  advantage  over  the  finite  difference 
technique . 

(4)  Continuum  derivatives 

Another  approach  for  determining  shape  design 
sensitivity  for  the  finite  element  method  stems  from  the  use 
of  the  concept  of  the  material  derivatives  in  conjunction 
with  the  adjoint  method.  The  approach  proceeds  by  forming 
material  derivatives  of  the  continuum  equation  [56,57]. 

Using  integration  by  parts  and  prescribed  boundary 
conditions,  the  derivatives  can  be  expressed  in  terms  of 
boundary  integrals,  which  is  a less  expensive  process  than  a 
direct  computation  of  the  derivative  dK/dx.  However,  when 
the  finite  element  method  is  used  for  analysis,  there  are 
many  numerical  difficulties  associated  with  the  evaluation 
of  boundary  integrals,  especially  for  low-order  curved 
boundaries . 

To  overcome  these  numerical  difficulties,  a domain 
method,  instead  of  a boundary  method,  is  developed  for  shape 
design  sensitivity  analysis  [53,58].  The  domain  integrals 
avoid  numerical  difficulties  associated  with  boundary 
integrals  but  they  are  as  expensive  to  use  as  the  implicit 
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differentiation  method  of  equation  (5.3). 

The  above  methods  are  all  well  developed  for  the  finite 
element  method.  With  the  exception  of  the  direct  difference 
approach  method  and  the  implicit  differentiation  method,  the 
other  three  methods  are  not  applicable  for  the  boundary 
element  method.  In  the  BEM,  the  matrices  [H]  and  [Q]  in 
equation  (3.20)  are  fully  populated  and  nonsymmetric , and 
efficient  solution  methods  that  are  applicable  to  banded  and 
symmetric  stiffness  matrices,  are  no  longer  applicable. 

More  recently,  Choi  and  Kwak  [32]  have  derived  a 
sensitivity  expression  involvinq  the  direct  boundary 
integral  equation  formulation  in  conjunction  with  the 
material  derivative  concept. 

In  the  following  section,  an  alternative  semi- 
analytical  approach  is  discussed  for  the  sensitivity 
analysis.  This  method,  which  utilizes  the  characteristics 
of  the  boundary  integral  equation  and  an  advantageous  use  of 
the  LU-decomposition  technique,  provides  an  efficient 
approach  for  the  task  at  hand. 

5 . 2 Semi-analytical  approach  in  BEM 

An  important  ingredient  in  shape  optimization  is  the 
computation  of  sensitivity  of  the  design  to  the  shape 
variables.  The  structural  shape  is  governed  by  a group  of 
independent  design  variables,  where  such  variables  affect 
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both  the  objective  function  and  the  structural  response  from 
which  the  design  constraints  are  computed. 

A careful  examination  of  the  eqn.(3.20)  indicates  that 
the  integral  functions  H and  U,  which  are  functions  of 
r= | £ - x|,  are  basically  determined  by  the  geometry  of  the 
structure.  When  one  of  the  shape  design  variables  is 
perturbed  by  a small  amount  6x,  only  a portion  of  the 
boundary  related  to  that  design  variable  will  be  changed.  It 
is  therefore  logical  to  expect  that  not  all  elements  in  the 
[H]  and  [G]  matrices  are  changed  during  a perturbation  of 
one  shape  variable.  This  concept  was  used  in  this  study  to 
improve  the  computational  efficiency  of  the  sensitivity 
analysis . 

A better  understanding  of  this  idea  is  obtained  by 
considering  the  boundary  of  a structural  domain  r, 
discretized  into  several  elements  as  shown  in  Figure  5.1. 

If  one  of  design  variables  d.  is  disturbed  to  d.+Sd.,  the 

1 l i' 

boundary  element  associated  with  the  design  variable  d^  will 
be  changed  to  r ' . From  the  integral  eqn.(3.16),  one  may 
note  that  if  the  load  point  x is  chosen  on  the  unchanged 
boundary,  the  difference  in  matrix  elements  of  [G]  and  [H] , 
and  [G']  and  [H'],  where  the  latter  represent  the  system 
matrices  of  the  perturbed  domain,  depends  on  the  location  of 
field  point  £ . That  is,  if  the  field  point  £ is  also 
located  on  the  unchanged  boundary,  these  matrix  elements 
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original 

boundary 


Figure  5.1  Original  and  perturbed  boundaries  discretized 
into  N segments. 
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will  not  be  changed.  Otherwise,  the  matrix  elements  related 
to  the  field  point  £ which  is  located  on  the  changed 
boundary,  must  be  recalculated.  Of  course,  if  the  load 
point  x lies  on  the  disturbed  boundary,  the  matrix  elements 
related  to  this  load  in  [H']  and  [G'J  will  also  change. 

For  simplicity,  one  may  consider  a linear  interpolation 
function  in  eguations  (3.18)  and  (3.19).  The  boundary  of 
the  structural  domain  is  discretized  into  N linear  segments 
as  shown  in  Figure  5.1.  The  system  eguations  of  the 
boundary  element  method  in  equation  (3.21)  can  be  written  as 
follows . 
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If  one  of  the  design  variables,  dj_,  is  perturbed  to 
d-L+Sd^,  boundary  elements  related  to  that  design  variable 
will  be  affected.  If  for  example,  the  design  variable 
affects  boundary  elements  between  those  denoted  by  indices 
K and  L-l  (see  Figure  5.1),  new  matrices  [A']  and  (C'}  are 
obtained  as  follows. 
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(5.5) 


In  the  equation  (5.5),  only  those  elements  within  the  cross- 
band are  changed  for  the  i-th  perturbed  design  variable. 

By  differentiating  equation  (3.21)  with  respect  to  the 
design  variable  d^,  the  following  relation  is  obtained, 

3x  3C  dA 

[A]  { } = { } - [ ] {x}  (5.6) 

ddi  ddi  3di 


where  the  matrix  [A]  on  the  left  hand  side  is  calculated 
before  the  design  variable  is  perturbed.  A finite  difference 
approximation  is  used  to  obtain  the  matrix  [dA/ddjJ  and  the 
vector  (ac/adi)  as  follows, 


3A 


[A']-[A] 


Ad, 


[ 


ad 


] 


(5.7) 
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dC  {C'}-{C> 

{ } = (5.8) 

3d  Adi 


where  [A']  and  (C'}  are  as  obtained  in  eqn.  (5.5),  requiring 
only  those  elements  of  [A]  and  { C } to  be  recalculated  that  are 
related  to  the  boundary  perturbations.  The  use  of  these 
finite  difference  approximations  is  the  reason  for  the 
approach  to  be  considered  as  a semi-analytical  approach. 

In  the  present  work,  the  system  equation  is  first  subjected 
to  an  LU-decomposition  as  follows. 

[A]  = [L]  [U]  (5.9) 

The  upper  and  lower  triangular  matrices,  [U]  and  [L] , 
respectively,  are  only  calculated  once  and  are  saved  for  all 
subsequent  design  variable  perturbations.  After  the 
sensitivities  [3A/3djJ  and  {SC/Sd^}  are  obtained,  the  solution 
(3x/3di)  of  eqn.  (5.6)  is  solved  in  two  stages.  First,  a 
forward  substitution  for  the  lower  triangular  matrix  of  [A] 
results  in  a solution  for  (H'}  as, 

[L]  (H'}  = (H)  (5.10) 

where , 


(H)  = (3C/3di)  - [ 3 A/  3 d jj  • (x) 

This  is  followed  by  a backward  substitution  for  the  upper 
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triangular  matrix  of  [A] . 

[U]  {ax/adi}  = { H ' } (5.11) 

In  the  approach  presented  above,  solution  of  the  system 
equation  for  each  perturbed  design  variable  to  obtain  the 
gradient  information,  is  replaced  by  the  backward  and  forward 
substitutions.  Since  the  matrix  [A]  is  fully  populated,  the 
efficiency  of  using  the  semi-analytical  approach  for 
sensitivity  analysis  is  evaluated  by  an  operation  count.  For 
a complete  set  of  gradient  evaluations,  the  comparison  of 
operation  counts  for  the  finite  difference  and  the  semi- 
analytical  approach  is  shown  in  Table  5.1.  In  comparing  with 
a direct  application  of  the  finite  difference  method,  the 
semi-analytical  approach  produces  computational  savings  of 
more  than  50%.  This  is  even  more  significant  in  the  presence 
of  a large  number  of  design  variables.  Additional  savings  in 
computational  cost  can  be  realized  by  selectively  computing 
those  elements  in  the  system  matrix  that  are  related  to 
perturbations  in  the  boundary. 
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Table  5.1  Comparison  of  operation  count  for  finite  difference 
and  semi-analytical  approaches  to  compute  design 
sensitivity. 


finite  difference  approach 

3 

N 2 

— ~ (NDV  + 1)  + N (NDV  + 1)  +0 (N) 

semi-analytical  approach 

3 

N 2 

— - + N (NDV  + 1)  + 0 (N) 

* [A]  is  an  N*N  matrix. 

* NDV  is  the  number  of  design  variable. 


CHAPTER  6 


GRID  REFINEMENT  AND  GRID  ADAPTATION 
IN  THE  BOUNDARY  ELEMENT  METHOD 


6 . 1 Introduction 


The  boundary  element  approach  offers  several  advantages 
over  the  finite  element  method  in  the  shape  design  problem. 
There  is,  however,  some  judgement  required  in  obtaining  an 
appropriate  discretization  of  the  boundary.  Since  the 
geometry  of  the  structural  domain  changes  continuously 
during  the  shape  design  process,  an  improper  grid 
distribution  would  yield  inaccuracies  in  stress  evaluation, 
particularly  when  the  geometry  is  complex. 

In  the  following  sections,  both  automated  grid 
refinement  and  grid  adaptation  methods  are  considered  to 
interface  with  the  optimum  shape  design  problem.  The 
concept  of  grid  refinement  technique  is  first  applied  to 
refine  the  grid  points  along  the  boundary  and  to  determine 
the  total  number  of  grid  points  for  the  structural  analysis. 
After  grid  refinement,  the  grid  points  are  redistributed 
along  the  structural  boundary  by  an  adaptive  scheme.  Two 


83 


84 


adaptive  schemes  are  studied  in  this  chapter.  The  first  is 
a variational  approach  which  provides  a simple,  explicit 
means  of  incorporating  desired  characteristics  into  the  grid 
generation  scheme.  The  second  approach  utilizes  existing 
optimization  algorithms  to  redistribute  the  grid  points  on 
the  structural  boundary. 


6 . 2 Grid  Refinement  Technique 

For  a two-dimensional  structural  domain,  the  grid 
distribution  on  the  boundary  is  usually  treated  as  a one- 
dimensional problem.  These  grid  points  not  only  determine 
the  nodes  at  which  analysis  is  performed,  but  also  define 
the  structural  boundary.  Hence,  proper  caution  must  be 
exercised  so  that  the  structural  geometry  is  not  altered 
when  redistributing  the  nodes  for  more  accurate  response 
analysis.  To  achieve  this  requirement,  a set  of  master 
nodes  are  introduced  along  the  structural  boundary.  The 
master  nodes  are  not  allowed  to  move  during  the  grid 
relocation  so  as  to  prevent  distortion  of  the  original 
domain. 

The  concept  of  grid  refinement  is  similar  to  the  h- 
method  of  FEM  mesh  refinement  [59].  The  grid  refinement 
considered  here  is  basically  a local  refinement  procedure. 
The  structural  boundary  is  assumed  to  consist  of  several 
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zones  which  are  defined  by  the  master  nodes.  On  the  basis 
of  a solution  computed  from  the  uniform  grid,  the  grid 
distribution  can  be  refined  locally  in  each  zone.  In  the 
proposed  approach  for  grid  refinement,  a control  function  is 
first  defined  and  computed  for  each  zone,  and  the  gradients 
of  this  control  function  determine  the  number  of  points  that 
must  be  removed  or  added  to  that  zone.  Once  the 
distribution  of  grid  nodes  is  defined,  the  total  number  of 
analytical  nodes  is  then  determined. 


6 . 3 Adaptive  Scheme  Based  on  Variational  Principles 

In  order  to  obtain  more  precise  analysis  results  in  the 
boundary  element  approach,  grid  points  need  to  be  clustered 
in  regions  of  large  response  gradients  [60].  In  this 
section,  a variational  approach  is  used  to  obtain  an 
adequate  grid  for  response  analysis.  The  desired  features 
of  such  a grid  are  represented  in  terms  of  a control 
function  defined  along  the  structural  boundary,  where  larger 
values  of  the  control  function  gradients  require  a 
proportionately  larger  number  of  computational  grid  points 
to  be  assigned  to  the  corresponding  region.  For  a problem 
pertaining  to  a one-dimensional  adaptation  of  grid  points 
along  the  boundary,  the  boundary  coordinate  s which  defines 
the  length  along  the  boundary,  and  a generalized  curvilinear 
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coordinate  £ which  pertains  to  the  computational  grid,  are 
as  shown  in  Figure  6.1.  One  form  of  a functional  Ip,  which 
can  be  used  as  a measure  of  adaptation  of  the  grid  spacing 
in  the  f direction,  can  be  written  in  terms  of  a control 
function  P as  [61], 

t 2 

s S 

ds  (6.1) 

P 

where  £s  = <3£/  ds,  and  is  referred  to  as  the  grid  density 
(number  of  grid  points  per  unit  length) . An  example  of 
choice  of  the  control  function  will  be  given  later.  For  the 
one-dimensional  problem,  the  relation  between  £s  and  s^  is 
expressed  as  follows. 


It  is  observed  from  eqn.  (6.1)  that,  for  the  functional  Ip 
to  be  minimum,  a region  on  the  boundary  with  a small  value 
of  control  function  P would  have  a corresponding  low  value 
of  £s.  Stated  differently,  a larger  grid  spacing  is 
obtained  in  the  region  of  small  P.  For  a boundary  with  a 
fixed  number  of  points,  the  above  implies  a denser  grid 
distribution  in  regions  of  high  P. 

The  variational  approach  yields  the  necessary 
conditions  for  optimality  that  a solution  for  the  grid 
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Figure  6.1  Definition  of  boundary  and  generalized 
curvilinear  coordinates. 
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density  must  satisfy  to  extremize  the  functional  Ip.  These 
are  the  familiar  Euler-Lagrange  equations,  and  for  the 
present  problem  result  in  a second  order,  one-dimensional 
Poisson  equation. 


?ss  - ^s  * Ps  / P = 0 (6'3> 

This  equation  governs  the  location  of  the  grid  coordinates 
in  the  physical  domain.  The  solution  of  this  equation  is 
facilitated  by  a transformation  to  the  curvilinear 
coordinate  system.  Using  a chain  rule,  the  second 
derivative  terms  can  be  shown  to  be  of  the  following  form. 


e 


ss 


(6.4) 


Substituting  equations  (6.2)  and  (6.4)  into  eqn.(6.3)  yields 
the  transformed  one-dimensional  equation  for  adaptation 
along  the  £ coordinate. 


+ * P£  7 P = ° (6'5) 

This  equation  can  be  solved  directly  for  the  grid  point 
spacing  along  the  £ -coordinate . Another  approach  of  looking 
at  the  adaptation  problem  is  obtained  if  one  integrates  eqn. 
(6.5)  once  to  obtain, 
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where, 


P • S£  = C 


(6.6) 


C 


I 


and,  represents  the  total  arc  length  along  the  boundary. 
It  is  obvious  from  the  above  expression  that  the 
relationship  between  the  control  function  P,  and  the  grid 
point  spacing  s , is  reciprocal. 

A proper  selection  of  the  control  function  P is  of 
essence  in  an  adaptive  relocation  of  grid  points  along  the 
structural  boundary.  In  the  problems  considered  in  this 
exercise,  the  control  function  was  defined  so  as  to  provide 
adequate  clustering  of  grid  points  in  regions  of  high  stress 
and  high  curvature  of  the  boundary.  One  such  control 
function  definition  is  of  the  following  form, 

f = £ * av  (1-£)  * ck  (6.7) 

where  f is  a control  function;  av  is  the  von  Mises  stress 
along  the  boundary;  c^  is  the  curvature  of  the  boundary;  p 
is  a weighting  factor  that  assumes  values  between  zero  and 
unity,  and  determines  the  relative  importance  of  the  stress 
or  curvature  in  the  adaptation  process.  The  ratio  of  the 
maximum  to  the  minimum  spacing  in  the  computational  grid  is 


90 


a measure  that  should  be  controllable  by  a proper  selection 
of  a control  function.  To  avoid  numerical  problems 
associated  with  the  control  function  assuming  a very  small 
value,  a modified  control  function  is  defined  as  follows, 

P = 1 + 7. f (6.8) 

where  the  function  f is  normalized  to  range  between  0 and  1, 
and  7 is  a smoothing  factor.  The  minimum  and  maximum  values 
of  P correspond  to  f=0  and  f=l,  respectively.  It  can  be 
shown  that  the  expression  for  grid  spacing  from  i to  i+1  is 
of  the  following  form. 


1 


t 1 + 7 • f i 

As  j_  = 


j 1 + 7 • f 

j 


(6.9) 


where  is  the  total  length  of  the  boundary,  and  P^  is  the 
value  of  the  control  function  between  i and  i+1,  and  can  be 
written  as  follows. 


Pi  - 1 + 7 . f i 


(6.10) 


Hence  the  ratio  of  the  maximum  to  the  minimum  spacing  may  be 
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expressed  as  follows. 


( ASi  ) 
( ASi  ) 


max 


min 


1 + 7 


(6.11) 


The  value  of  this  ratio  is  clearly  dependent  on  the 
smoothing  factor  7 which  is  prescribed  in  the  control 
function.  However,  as  explained  in  subsequent  sections,  the 
approach  used  in  this  work  for  grid  adaptation  and 
refinement  is  relatively  insensitive  to  the  choice  of  this 
parameter. 


6 . 4 Grid  Adaptation  Obtained  by 

Nonlinear  Programming  Methods 

The  variational  approach  described  in  the  previous 
section  to  obtain  the  computational  grid  is  applicable  in 
several  engineering  problems.  However,  if  the  functional  Ip 
in  equation  (6.1)  is  complicated,  or  if  additional 
constraints  must  be  imposed  in  the  grid  adaptation,  a 
solution  for  grid  point  spacing  is  not  as  easily  available. 

The  use  of  nonlinear  programming  methods  to  relocate  the 
grid  points  is  a more  generally  applicable  strategy.  This 
approach  can  be  extended  to  include  additional  constraint 
conditions.  Furthermore,  in  the  event  that  several 
noncommensurate  criteria  are  to  be  considered  for  the  grid 
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point  locations,  a multicriterion  optimization  approach  can 
be  adopted.  Instead  of  solving  the  partial  differential 
equation  for  the  optimal  grid  point  locations,  a nonlinear 
programming  algorithm  is  utilized  to  minimize  the  functional 
Ip.  If  the  boundary  of  the  structure  is  discretized  into  N 
segments,  the  functional  Ip  can  be  written  as  follows, 

N 1 

Ip  = S (6.12) 

i=l  Pi  Asi 

where,  as  before,  Pi  is  the  control  function  in  the  i-th 
region  and  Asi  represents  the  length  between  two  adjoining 
nodes,  (i-1)  and  (i).  The  mathematical  statement  of  the 
nonlinear  programming  based  optimization  problem  for  grid 
adaptation  can  be  written  as  follows. 

Find  Si,  S2  , ...  , sjj 

to  minimize  Ip  (6.13) 

Here  s^,s2,  ...  , s^  are  the  length  coordinates  of  all  moving 

grid  points.  It  is  important  to  note  that  this  set  does  not 
include  the  master  nodes.  This  statement  is  basically  that 
of  an  unconstrained  optimization  problem,  for  which  several 
efficient  solution  algorithms  such  as  the  Fletcher-Reeves 
method  [62]  and  the  Broydon-Fletcher-Goldf arb-Shanno 
variable  metric  method  [63-66]  are  available.  The  gradient 
of  the  objective  function  is  needed  in  such  first  order 
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methods,  and  can  be  easily  obtained  by  forward  or  backward 
finite  difference  method.  However,  the  size  of  the 
perturbation  step  in  such  a difference  approach  is  very 
important.  Since  the  master  nodes  are  not  allowed  to  move 
during  the  grid  redistribution,  an  arbitrary  perturbation  of 
the  design  variables  can  cause  a grid  point  to  be  located  in 
an  undesired  zone.  This  is  most  likely  to  happen  if  the 
perturbation  is  a fixed  fraction  of  the  design  variable,  as 
is  frequently  the  case  in  member  sizing  problems.  In  the 
present  work,  a fixed  design  variable  perturbation  was 
applied  to  all  design  variables. 

The  problem  of  grid  adaptation  for  more  than  one 
criterion  function,  and  subsequence  use  of  the  nonlinear 
programming  approach,  is  explained  as  follows.  Consider  a 
control  function  that  can  be  separated  into  two  distinct 
objective  criteria,  one  based  on  the  behavior  response  such 
as  the  von  Mises  stress,  and  a second  based  on  the  local 
curvature  of  the  structural  boundary.  The  location  of  grid 
points  must  be  optimal  from  the  standpoint  of  both  criteria, 
resulting  in  a multicriterion  optimization  problem.  In  this 
study,  an  implementation  of  the  global  criterion  approach 
[67]  was  used  to  obtain  solutions  to  the  grid  adaptation 
problem.  This  approach  is  summarized  here  for  completeness. 

Consider  a situation  in  which  the  control  function 
consists  of  two  candidate  criteria,  fi(s)  and  f2(s), 
representing  the  importance  of  the  von  Mises  stress  function 
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(ctv)  and  the  curvature  of  the  boundary  (ck) , in  the  problem 
of  locating  grid  points.  These  functions  can  be  written  as 
follows . 


N 1 

f 1 (s)  = s 

i=l  aV i ASi 


(6.14) 


N 1 

f2(s)  = S 

i=l  ck-[  Asi 


(6.15) 


Each  candidate  criterion  may  be  optimized  individually, 

ignoring  the  existence  of  the  second  criterion.  The 

resulting  optimal  values  of  the  objective  functions  filc^(s) 
i d 

and  f2  (s)  are  referred  to  as  the  ideal  solutions  for  each 
criterion  function.  The  global  criterion  method  is  based  on 
the  premise  that  the  design  vector  for  which  both  criteria 
are  simultaneously  optimized  may  be  obtained  by  minimizing  a 
metric  measure  of  the  distance  between  the  ideal  solutions 
and  the  true  optimum.  As  shown  in  [67],  this  results  in  a 
min-max  problem  which  can  be  converted  into  an  equivalent 
scalar  minimization  problem  as  follows. 

Minimize  p 
Subject  to 

f ! (s)  - fild(s) 


- p < 0 
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f2(s)  ~ f2ld(s) 

w2  Id 

f2  (s) 


- p < 0 


where , 


w 1 + w 2 = 1 


Here  «]_  and  w 2 are  the  weighting  coefficients  representing 
the  relative  importance  of  each  criterion.  Figure  6.2 
illustrates  the  feasible  and  infeasible  spaces  for  a two- 
criterion  function.  By  choosing  different  weighting  factors 
and  ^2,  the  ideal  solution  fld(s)  is  located  on  the 
boundary  between  f1ld(s)  and  f2ld(s) . 


6 . 5 Examples  of  Grid  Refinement  and  Adaptation 


The  primary  motivation  for  using  an  adaptive  grid  is  to 
allocate  an  appropriate  number  of  grid  points  to  regions  of 
high  stress  and  curvature.  In  problems  of  shape  design 
using  the  BEM,  grid  points  play  a dual  role  of  defining  the 
structural  shape  and  providing  discrete  points  for  analysis. 
In  computing  the  sensitivity  of  the  structural  response  to 
changes  in  the  boundary,  where  the  latter  is  linked  to  the 
position  of  grid  points,  proper  care  must  be  exercised  to 
eliminate  distortions  in  the  geometry  of  the  boundary.  The 
master  nodes  uniquely  define  the  shape  of  the  structure  and 
serve  as  points  for  response  analysis  as  well  as  provide  the 
interpolation  points  for  the  structural  geometry  definition. 
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Figure  6.2  Graphical  representation  of  a two-criterion 
function  design  space. 
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These  nodes,  shown  in  Figure  6.3(a),  are  not  moved  when  the 
adaptive  scheme  is  applied. 

Grid  adaptation  and  refinement  is  essentially  a two- 
stage  process.  A uniformly  distributed  grid  is  first 
assumed  between  the  master  nodes,  and  the  corresponding 
approximate  solution  is  determined.  This  solution  allows  a 
computation  of  the  control  function  and  its  gradient  along 
the  boundary.  The  approach  adopted  is  one  in  which  more 
grid  points  are  assigned  to  the  high  gradient  regions,  and 
less  or  no  grid  points  are  located  in  regions  with  a smooth 
variation  of  the  control  function.  A relatively  simplistic 
scheme  was  implemented  to  realize  this  goal.  The  numerical 
values  of  the  gradients  of  the  control  function  were 
classified  into  several  levels,  and  a proportional  number  of 
grid  points  were  assigned  to  each  of  these  levels.  The 
effect  of  such  a redistribution  of  points  based  on  control 
function  gradients  is  illustrated  in  Figure  6.3(b). 

Once  the  total  number  of  grid  points  for  each  region  is 
determined,  the  governing  equation  for  grid  adaptation  is 
employed  to  provide  a better  grid  redistribution.  This  is 
shown  in  Figure  6.3(c).  Note  that  in  the  solution  of  grid 
adaptation,  the  master  nodes  are  always  fixed  to  maintain 
the  boundary  of  the  structure. 

Figure  6.4  illustrates  the  process  of  grid  point 
distribution  on  the  basis  of  control  function  values.  After 
an  approximate  solution  is  obtained  for  a uniform  grid 
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distribution  as  shown  in  Figure  6.4(a),  the  grid  is  refined 
between  each  pair  of  neighboring  master  nodes  on  the  basis 
of  the  gradient  of  the  control  function.  In  Figure  6.4(a), 
we  observe  that  the  gradients  of  the  control  function  around 
points  A,  B and  C are  much  higher  than  elsewhere.  Thus, 
several  new  grid  points  are  assigned  to  these  regions  and 
evenly  distributed  within  the  regions.  A smaller  number  of 
grid  points  are  assigned  to  regions  of  low  gradients,  as 
shown  in  Figure  6.4(b).  The  points  placed  in  these  regions 
are  next  relocated  by  a solution  of  the  adaptive  equation, 
and  the  results  of  this  adaptation  are  shown  in  Figure 
6.4(c).  It  is  worthwhile  to  bear  in  mind  that  since  master 
nodes  have  been  introduced  along  the  boundary,  the  maximum 
spacing  (^si)max  •*-n  the  eqn.  (6.11)  is  equal  to  the  space 
between  master  nodes.  This  also  renders  the  adaptation 
scheme  relatively  insensitive  to  the  choice  of  the  smoothing 
factor  7 . 

Table  6.1  compares  the  von  Mises  stress  distribution 
along  a structural  boundary  for  a problem  solved  with  149 
evenly  distributed  analytical  points,  and  with  61  analytical 
points  assigned  on  the  basis  of  grid  refinement  and 
adaptation.  The  differences  in  the  numerical  results  are 
small,  even  though  the  number  of  grid  points  in  the  second 
case  are  less  than  half  the  number  in  the  first.  Such  a 
reduction  of  grid  points  can  result  in  substantial  savings 
in  computational  cost  in  shape  design  problems. 


bSSyfUnCtl°n  dist«bution  along  the 

■4<b>  control ffunction^aSe<^  °n  gradients  the 


Control  Function 
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Figure  6.4(c)  Grid  adaptation  in  regions  of  high  stress. 
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Table  6.1  Comparison  of  ctv/(7  ,,  between  149  uniform  grids  and  61 
redistributed  grids. 


selected 
master  nodes 

av/aal] 

case  1 : 

uniform  grid  with 
149  analytical  nodes 

case  2 : 

adaptive  grid  with 
61  analytical  nodes 

1 

0.8500 

0.8629 

2 

0.1207 

0.1197 

3 

0.0999 

0.1029 

4 

0.1144 

0.1121 

5 

0.1234 

0.1185 

6 

0.1144 

0.1121 

20 

0.5154 

0.5092 

21 

0.0596 

0.0595 

22 

0.0493 

0.0538 

23 

0.0650 

0.0480 

24 

0.0411 

0.0404 

25 

0.0259 

0.0254 

26 

0.0019 

0.0096 

27 

0.0470 

0.0484 

28 

0.1166 

0.1144 

29 

0.2640 

0.2690 

* a v is  von  Mises  stress. 

* CTan  is  allowable  stress  (62,000  N/cm2  ) . 
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Another  example  illustrating  grid  adaptation,  based  on 
a nonlinear  programming  solution  of  a multicriterion 
optimization  problem,  is  presented  next.  Hypothetical 
variations  for  av  and  ck,  defined  in  preceding  sections, 
were  assumed.  These  can  be  written  as  explicit  functions  of 
the  coordinate  variables  as  follows. 

av(s)=sin(s)-0.6*sin(s)+0.1*sin(s) 

ck(s)=  |sin(3*s) | 

Figure  6.5  shows  a distribution  of  nineteen  nodes  in  the 
interval  between  0 and  n , with  the  two  end  points  fixed. 
Figures  6.5(a)  and  6.5(b)  show  the  adapted  grid 
distributions  based  on  av(s)  and  ck(s) , respectively.  The 
values  of  the  ideal  solutions  f^ld(s)  and  f2ld(s)  were 
113.44  and  109.09,  respectively.  Figure  6.5(c)  shows  the 
grid  distribution  when  crv(s)  and  ck(s)  were  considered 
egually  important  in  the  adaptation  process.  This 
distribution  was  obtained  as  a solution  of  a multicriterion 
problem  where  the  weighting  factor  u>±  and  o>2  were  both  set 
egual  to  0.5.  Figure  6.5(d)  shows  the  grid  distribution  for 
a case  where  curvature  is  weighted  more  than  the  stress 
distribution.  This  distribution  corresponds  to  values  of 
ui=0. 333  and  ^2=0. 667. 
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Figure  6.5  Grid  distribution  obtained  by  using  a control 
function  of  the  form  + u2-ck. 


case 

case 

case 

case 


(a) 

(b) 

(c) 

(d) 


wi  = w2  = o, 

“1  = 0,  w2  = 1, 

w]_  = 0.5,  w2  = 0.5, 

— 0.333,  w2=0.667, 
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6 . 6 Optimization  examples  using  grid  refinement 
and  grid  adaptation  techniques 

In  this  section,  two  optimization  examples  are  used  to 
illustrate  the  function  and  the  efficiency  of  applying  the 
grid  refinement  and  grid  adaptation  schemes  in  conjunction 
with  BEM  for  shape  optimal  design.  The  framework  of  the 
optimization  procedure  is  shown  in  Figure  6.6.  The  semi- 
analytical  approach  is  employed  to  obtain  the  gradient 
information . 

A flat  plate,  supported  as  shown  in  Figure  6.7,  is 
subjected  to  a concentrated  load  of  45,000  N.  The  allowable 
stress,  Poisson  ratio,  and  Young's  modulus  for  the  plate 
material  are  620  Mpa,  0.3,  and  226  Gpa,  respectively.  The 
dimensions  of  the  plate  and  a definition  of  the  shape  design 
variables  is  also  shown  in  this  figure.  The  distance  between 
the  support  is  fixed,  as  are  the  points  A and  B.  Other 
boundaries  are  defined  in  terms  of  half-parabolic  curves. 
Design  variables  x^,  X2  and  X3 , X4  control  the  right  outer 
and  the  left  outer  boundaries,  respectively.  Design 
variables  X5  and  xg  define  the  height  of  the  two  parabolic 
curves  of  the  inner  boundary.  Symmetry  of  the  loading  and 
support  dictate  that  x^=X3  and  X2=X4 . 

Figure  6.8  illustrates  the  design  iteration  history  for 
this  plate.  Starting  from  the  initial  design  of  the 
structure  shown  in  Figure  6.8(a),  the  design  after  four 
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System  equations  are  solved  by 
LU-decomposition . 

Upper  and  lower  triangular  matrices 
are  saved  for  sensitivity  evaluation. 
Compute  objective  function  and  constraint 
information. 


Use  the  semi-analytical  approach 
to  evaluate  the  gradient  of  the 
objective  function  and  constraints. 
Forward  and  backward  substitution 
is  used  to  solve  the  system  equations. 


Figure  6.6  A fully  nonlinear  optimization  algorithm  with  grid 
ridistribution  and  semi-analytical  computation  of 
design  sensitivity. 
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45  KN 


Figure  6.7  Definition  of  geometry  and  design  variables  for 
the  elastic  plate  problem. 
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cycles  is  as  shown  in  Figure  6.8(b).  At  this  point,  the 
original  node  distribution  is  no  longer  suitable  for 
analysis,  and  the  application  of  grid  refinement  and  grid 
adaptation  results  in  a node  distribution  shown  in  Figure 
6.8(c).  The  optimum  shape  of  the  structure  is  obtained  in 
the  sixth  cycle,  and  is  shown  in  Figure  6.8(d). 

The  second  test  problem  involved  the  minimum  weight 
design  of  a fillet  to  sustain  a uniaxial  distributed  load  of 
40,000  N/cm,  as  indicated  in  Figure  6.9.  The  material 
properties  are  the  same  as  in  the  previous  example.  The 
structural  boundary  to  be  redesigned  is  included  between 
points  A and  B,  as  shown  in  the  Figure  6.9.  This  boundary 
was  discretized  into  ten  equal  segments,  and  all  nodes  on 
this  line  (excluding  points  A and  B)  were  allowed  to  move 
perpendicular  to  the  line  AB  during  the  redesign  process. 
This  results  in  a total  of  nine  design  variables  for  the 
problem. 

This  is  a classical  problem  in  optimal  shape  design,  and 
is  known  to  present  problems  in  obtaining  a converged 
optimum.  Changes  in  boundary  curvature,  particularly  the 
introduction  of  sharp  corners  on  the  boundary,  results  in 
significant  movements  of  the  zones  of  high  stress.  The 
numerical  optimization  exploits  inherent  weaknesses  in 
discrete  analysis  methods  related  to  the  fact  that  analysis 
is  only  performed  at  predetermined  points.  The  design 
process,  therefore,  exhibits  an  oscillatory  behavior.  There 


109 


■p 

(0 

43 

-p 

QJ 

• c 

p 

•P 

g 0 

3 

3 

<U  *H 

0 

iH 

rH  -p 

<p 

ft 

43  3 

0 13 

>p 

0 

P "P 

0 

•P 

ft  P 

■P 

•P 

Cn 

w 

<u  w 

c 

3 

•P  -P 

•P 

ip 

(0  T3 

c 

(1) 

P • 

•P 

ftT3  a) 

Cn 

a) 

•P  i — 1 

Q) 

13 

0 P 0 43 

■p 

•h  trN 

-P  0 -P 

P 

[/)  iH 

<v 

0 

3 (0  13 

<w 

rH  -p  -p 

G 

<U  -P  P 

0 

c 

•H  3 

•H 

tr> 

P G 0 -P 

•p 

0 -P  <W 

3 

V) 

<p 

43 

a) 

C <D 

•H 

T3 

c o £ 

p 

cn  -p  .p 

1 — I 

-H  T} 

W 

3 

w <d  <p 

•H 

E 

(D  M Ofl 

■P 

•a  3 

-P 

13  CPTS 

ft  • 

'-I  G • 

■H 

0 E 

3 G -P 

U • 

<u 

•H  Cr>  G 

Cr>  (U  r-H  i— i 

■P  -H  -P 

rH 

3 13 

•h  tn  tr>  5 u 

G 0 

G <D  <1) 

0)  Jj 

HQ13! 

2 O 

ft  ft 

<0  43 

O 

T? 



CO  CO  CO 

CO 

VO  VO  KO 

• 

VO 

<D 

P 

3 

•P 

ft 


110 


E 

u 


Figure  6.9  Definition  of  geometry  and  loading  for  the 
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are  essentially  two  approaches  to  counter  oscillations  in 
such  problem.  The  first  espouses  the  use  of  a very  large 
number  of  points  on  the  boundary  to  obtain  better  resolution 
of  regions  of  high  stress.  This  has  the  associated  drawback 
of  larger  computational  requirements.  The  other  approach  is 
one  which  uses  grid  refinement  and  adaptation  techniques  to 
locate  a fewer  number  of  points  in  the  critical  regions. 

The  order  of  the  resulting  system  equations  is  low,  yielding 
superior  computational  efficiency. 

The  initial  design  for  the  fillet  problem  is  as  shown  in 
Figure  6.10(a).  An  adaptive  grid  was  defined  and  the  design 
allowed  to  proceed  with  prescribed  move  limits  on  each 
design  variable  during  the  optimization.  The  presence  of 
such  move  limits  alleviates  the  computational  burden  to  some 
extent,  as  a new  grid  refinement  is  not  necessary  after  each 
iteration.  Figures  6.10(b)  and  6.10(c)  illustrate  the 
progression  of  the  design  to  an  optimum. 
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Initial  design  for  the  fillet  problem. 
Design  at  begining  of  third  cycle. 

Final  optimal  design  for  the  fillet  problem 
obtained  at  the  fifth  cycle. 


CHAPTER  7 


AN  INTEGRATED  APPROACH  FOR  SHAPE  OPTIMIZATION 
USING  THE  BOUNDARY  ELEMENT  METHOD 


7 . 1 Introduction 


The  use  of  an  integrated  formulation  in  structural 
optimization  was  first  proposed  by  Schmit  and  Fox  [5,6]. 

The  method  proposed  collapsing  all  constraints,  including 
the  equality  constraints  (a  set  of  linear  equilibrium 
equations)  and  inequality  constraints  (behavior 
constraints) , into  an  unconstrained  function  of  several 
design  variables,  and  obtaining  a solution  to  the 
unconstrained  problem  by  the  penalty  function  method.  This 
approach  was  abandoned  at  the  time  for  lack  of  efficient 
methods  to  solve  the  poorly  conditioned  unconstrained 
optimization  problem. 

In  view  of  the  recent  developments  in  nonlinear- 
programming-based  optimization  methodology,  the  approach  of 
simultaneous  analysis  and  design  has  again  surfaced  as  a 
viable  option.  Haftka  [7],  uses  an  element-by-element 
preconditioned  conjugate  gradient  method  in  the  solution  of 
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a linear  analysis  problem,  and  Newton's  method  in  a design 
subject  to  a nonlinear  collapse  load  constraint. 

When  optimization  techniques  are  integrated  with 
structural  analysis,  the  design  problem  can  be  viewed  as  a 
nested  optimization  problem.  Since  the  computation  of 
behavior  variables  such  as  stresses  and  displacements  is 
unnecessary  for  each  trial  design  vector  x,  the  use  of  this 
approach  reduces  the  solution  of  the  analysis  equations  to 
an  iterative  process,  which  proceeds  in  conjunction  with  the 
optimization  solution.  This  approach  not  only  eliminates 
the  evaluation  of  a large  number  of  intermediate  trial 
designs,  but  also  achieves  optimum  results  of  both  the 
response  solution  and  shape  design  simultaneously. 

Most  integrated  approach  formulations  available  in  the 
literature  pertain  to  the  structural  member  sizing  problems, 
and  are  limited  to  finite  element  based  analysis.  In  this 
study,  the  integrated  approach  is  extended  to  shape 
optimization  problems,  wherein  the  BEM  analysis  equations 
are  represented  as  a set  of  equality  constraints  in  the 
optimization  problem.  This  approach  precludes  an  explicit 
solution  of  the  set  of  linear  algebraic  equations  (BEM 
equations) , at  the  cost  of  an  increased  number  of 
constraints  and  design  variables  (response  variables  are  now 
design  variables)  in  the  optimization  problem.  The 
generalized  reduced  gradient  (GRG)  method  is  employed  to 
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obtain  a solution  to  the  combined  analysis  and  design 
problem. 

7 • 2 Integrated  Approach  Using  The  Boundary  Element  Method 

In  the  conventional  structural  optimization  approach, 
the  optimization  and  analysis  techniques  are  independent 
tools  in  the  design  process.  The  optimization  approach 
provides  a mathematical  tool  to  effectively  search  the 
design  space  so  as  to  achieve  an  optimal  objective  function 
value,  and  simultaneously  satisfy  a set  of  constraints.  The 
analysis  simply  provides  the  optimizer  with  necessary 
information  to  proceed  towards  the  optimum.  In  this  study, 
an  integrated  formulation  is  adopted,  coupling  the  boundary 
element  analysis  equations  with  the  nonlinear  programming 
based  optimization  problem. 

The  BEM  formulation  for  analysis  of  a structure 
typically  leads  to  a system  of  linear  equations, 

[A]  {v}  = {C}  (7.1) 

where  the  response  quantities  v in  equation  (7.1),  are 
generally  solved  by  Gaussian  elimination.  In  the  integrated 
approach,  these  linear  equations  are  treated  as  additional 
equality  constraints  in  the  optimization  problem,  combining 
the  components  of  the  unknown  response  vector  v with  the 
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vector  of  shape  design  variables  d.  Further,  the  equations 
relating  the  stress  response  and  the  response  vector  v are 
also  included  as  additional  equality  constraints.  Such 
equality  constraints  relate  the  von  Mises  stress  to  the 
response  variables  v and  are  represented  as, 

ai  ~ CTVi(v)=  0 (7.2) 

where  cr-^  is  the  stress  design  variable  which  represents  the 
value  of  von  Mises  stress  at  boundary  node  i,  and  avi(v) 
represents  the  von  Mises  stress  as  a function  of  the 
response  variables  v at  node  i.  Additional  inequality 
constraints  are  obtained  due  to  requirement  of  allowable 
limits  on  the  stress  response. 


CTi 


"iall  s 0 


(7.3) 


where  CTial^  is  the  i-th  allowable  stress  prescribed  at  the 
beginning  of  the  design  cycle.  A general  mathematical 
statement  for  the  integrated  optimization  problem  involving 
the  boundary  element  method  may  be  written  as  follows. 
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Find  d,  v,  and  to  minimize  W 


Subject  to  Av  - C = 0 


a j_  - crvi(v)  = 0 


° i 


aiall 


(7.4) 


where  d is  a vector  of  shape  design  variables,  and  W is  the 
structural  weight.  For  a given  set  of  design  variables,  the 
structural  weight  is  only  a function  of  shape  design 
variables  d.  The  matrices  A and  C in  the  equality 
constraints,  and  the  boundary  stresses  cri(v)  in  the 
inequality  constraints  are  obtained  from  the  boundary 
element  formulation. 

As  was  stated  in  Chapter  4,  numerical  accuracy  of  the 
stress  response  is  extremely  sensitive  to  the  grid 
distribution.  Improper  node  point  locations  would  result  in 
misleading  information  from  the  BEM  based  analysis. 

Although,  a uniformly  dense  grid  distribution  for  the  entire 
boundary  may  yield  better  accuracy,  the  process  can  become 
computationally  prohibitive.  In  using  an  integrated 
optimization  strategy  for  a two-dimensional  elasticity 
problem,  the  addition  of  one  node  results  in  the 
introduction  of  three  equality  constraints  (two  for  boundary 
element  equations  and  one  for  stress  equality  equation)  and 
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one  additional  inequality  constraint.  Furthermore,  three 
more  design  variables  (two  v and  one  ) are  added  into  the 
design  space.  From  the  standpoint  of  computational 
efficiency  and  reduction  of  memory  requirements,  the  total 
number  of  nodes  must  be  minimized.  The  proper  node 
distribution  required  in  the  solution  of  the  problem  is 
determined  by  the  grid  refinement  and  adaptation  techniques 
discussed  in  the  Chapter  6. 

The  shape  design  variable  vector  d changes  continuously 
during  the  optimization  process.  After  just  a few 
iterations,  the  initial  grid  distribution  may  become 
unsuitable  for  the  new  boundary  shape.  To  avoid  a new  grid 
definition  for  every  iteration,  move  limits  were  imposed  on 
the  shape  design  variables.  A prescribed  grid  distribution 
can  therefore  remain  effective  for  an  admissible  and 
controlled  variation  of  the  boundary. 

The  adaptive  grid  distribution  method  was  used  in 
conjunction  with  an  integrated  analysis/design  procedure  to 
determine  the  minimum  weight  geometry  of  a fillet  to  sustain 
a uniaxial  distributed  load  of  35,000  N/cm.  This 
minimization  was  subject  to  constraints  of  maximum  allowable 
value  of  the  von  Mises  equivalent  stress  on  the  boundary. 

The  dimensions  and  material  properties  of  the  structure  are 
indicated  in  Figure  7.1.  As  shown  in  this  figure,  the 
structural  boundary  to  be  redesigned  is  included  between 
points  A and  B.  This  boundary  was  discretized  into  eight 
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Figure  7.1  A planar  fillet  design  problem. 
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equal  segments,  and  all  nodes  on  this  line  (excluding  points 
A and  B)  were  allowed  to  move  perpendicular  to  the  line  AB 
during  the  redesign  process.  This  results  in  a total  of 
seven  design  variables  for  the  problem.  After  grid 
refinement  and  adaptive  procedures,  total  number  of  grid 
points  and  the  grid  locations  are  determined  for  structural 
analysis.  In  this  problem  a total  of  37  nodes  were  retained 
for  structural  analysis.  These  nodes  were  distributed  as 
shown  in  figure  7.2.  The  number  of  design  variables  for  the 
integrated  approach  is  121,  which  consists  of  7 shape  design 
variables  d,  74  response  design  variables  v,  and  37  stress 
design  variables  a 

The  values  of  components  of  the  vector  v for  this 
problem  typically  ranged  from  10  3 to  10  7 . An  improper 
choice  of  initial  values  of  v resulted  in  inordinately  high 
computational  expense,  and  may,  in  some  cases,  preclude  the 
determination  of  an  optimum  design.  A rational  choice  of 
initial  design  variable  values  is  therefore  essential  in  the 
integrated  approach.  For  a given  set  of  initial  shape 
variables  d,  exact  value  estimates  of  the  response  design 
variable  vector  v and  aj  can  be  obtained  by  an  evaluation 
the  BEM  based  analysis  equations.  This  process  can  be 
periodically  repeated  in  the  iterative  design  process. 

The  total  number  of  constraints  is  dependent  on  the  number 
of  grid  points.  In  the  present  problem,  there  were  74 
equality  constraints  corresponding  to  the  linear  equations 
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Figure  7.2  37  adaptatively  redistributed  node  points 

for  fillet  problem. 
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(7.1),  37  equality  constraints  for  the  stress  response 
computations  (7.2),  and  37  inequality  constraints  for  stress 
allowable  limits  (7.3),  respectively.  Since  the  boundary 
element  equations  are  linear  with  respect  to  the  response 
variable  vector  v,  the  matrix  A can  be  directly  used  as  the 
gradient  of  the  equality  constraints  corresponding  to  the 
BEM  analysis  equations.  The  gradient  information  for  all 
other  constraints  may  be  obtained  by  the  finite  difference 
approach. 

The  integrated  optimization  problem  was  solved  by  the 
generalized  reduced  gradient  (GRG)  method,  with  prescribed 
move  limits  established  for  each  shape  design  variable.  The 
changes  of  shape  design  variables  for  each  piecewise  linear 
cycle  were  restricted  to  within  20%  of  starting  values. 
Within  these  move  limits,  first  order  linear  approximations 
of  the  constraints  were  used.  The  corresponding  results 
were  obtained  in  five  cycles  and  the  final  shape  is  shown  in 
Figure  7.3.  Such  an  integrated  approach  eliminates  the 
necessity  for  repeated  decompositions  of  the  boundary 
element  equations,  as  would  be  necessary  in  the  traditional 
formulation.  However,  it  is  worthwhile  to  note  that  the 
integrated  approach  results  in  an  increase  in  design 
variables  and  constraints,  giving  rise  to  an  increase  in 
computational  costs.  This  increase  is  largely  due  to  the 
fact  that  tight  move  limits  must  be  imposed  in  the 
optimization  problem  because  of  the  poor  numerical  scaling 
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Figure  7.3  Final  shape  design  for  fillet  problem. 
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of  the  design  space.  For  instance,  the  value  of  equality 
constraints  (7.1)  must  necessarily  be  very  small  to  obtain 
accurate  values  of  the  response  variables.  Multiplying 
factors  of  order  102  to  103 , were  used  to  scale  the  equality 
constraints  to  improve  the  accuracy  of  the  solution  of  the 
response  design  variables,  stress  response  quantities  (a^) , 
and  the  shape  design  variables  (d) , in  the  integrated 
approach. 


CHAPTER  VIII 


CONCLUSIONS  AND  RECOMMENDATIONS  FOR 
FUTURE  RESEARCH 


Some  key  ideas  related  to  an  efficient  optimal  shape 
synthesis  of  elastic  structures  have  been  discussed  in  this 
study.  The  boundary  element  method  is  used  for  analysis  as 
it  provides  better  accuracy  in  response  computations  at  the 
boundary.  Two-dimensional  elastic,  torsional  shaft 
problems,  and  small-strain  elastostatics  problems  are 
successfully  solved  by  the  BEM.  Both  constant  and  linear 
elements  are  used  for  the  boundary  element  analysis 
necessary  for  shape  synthesis. 

A strategy  for  choosing  proper  shape  design  variables  is 
examined.  The  method  is  based  on  the  computed  load  path  in 
the  structural  domain.  Gradient-based,  nonlinear 
programming  algorithms  are  investigated  in  the  shape 
optimization  problems.  In  general,  the  direct  NLP  methods 
require  fewer  computations  of  objective  function  and 
gradients  than  the  indirect  NLP  method.  For  the  NLP  problem 
involving  several  equality  constraints,  such  as  in  the 
integrated  approach,  the  GRG  method  is  effectively  employed 
to  obtain  optimum  solutions. 
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The  evaluation  of  gradients  in  the  NLP  method  is  done  by 
semi-analytical  methods,  where  such  methods  are  shown  to  be 
significantly  more  efficient  than  finite  difference 
approximations.  In  the  present  work,  this  approach  was 
implemented  in  computation  of  design  sensitivity  for  a 
torsional  shaft  problem  and  for  a class  of  two  dimensional, 
linear  elastostatics  problem.  The  concept  of  this  gradient 
evaluation  can  be  readily  extended  to  high  order  elements  of 
BEM  for  shape  optimization  problems. 

The  research  also  describes  an  effective  approach  for 
grid  refinement  and  adaptation  in  shape  synthesis  problems. 
The  minimization  of  the  total  number  of  grid  points  and 
their  proper  location  is  especially  significant  in  the  BEM 
formulation.  This  is  true  from  the  standpoint  of  both 
numerical  accuracy  and  computational  efficiency.  The  total 
number  of  analytical  nodes  is  determined  by  grid  refinement. 
Two  adaptive  schemes,  one  based  on  a variational  approach, 
and  the  other  on  a nonlinear  programming-based  unconstrained 
optimization  algorithm,  are  examined  in  this  work.  It  is 
worthwhile  to  note  that  the  nonlinear  programming  method 
allows  for  addition  of  constraints  as  well  as  an  ability  to 
account  for  multiple  criteria  in  the  grid  adaptation.  In 
using  such  grid  refinement  and  adaptive  schemes,  the  order 
of  system  equations  in  the  BEM  analysis  can  be  reduced 
significantly.  To  avoid  a new  grid  definition  for  every 
iteration,  move  limits  were  imposed  on  the  shape  design 
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variables.  A prescribed  grid  distribution  can  therefore 
remain  effective  for  an  admissible  and  controlled  variation 
of  the  boundary. 

The  research  also  examines  the  feasibility  of  using  an 
integrated  formulation  for  a BEM  based  optimum  shape 
synthesis.  The  efficiency  of  the  proposed  scheme  is 
attributed  to  eliminating  the  necessity  of  decomposing  the 
system  matrix  in  the  BEM  formulation.  This  matrix  is  not 
banded  as  in  the  finite  element  formulation,  and  is 
therefore  not  amenable  to  efficient  decomposition  methods. 
The  inclusion  of  behavior  variables  as  additional  design 
variables  in  the  integrated  optimization  formulation  often 
results  in  a poorly  scaled  system,  and  detracts  from  the 
efficiency  of  the  process. 

With  regard  to  future  work  in  this  direction,  a 
substructure  shape  optimization  concept  which  combines  both 
the  FEM  and  the  BEM  is  highly  recommended.  In  many 
practical  design  problems  encountered  in  automotive  or 
aircraft  structures,  the  structural  domain  is  very  large, 
and  can  not  be  properly  discretized  by  the  boundary  element 
method.  The  FEM  can  be  first  applied  to  analyze  the  entire 
structural  domain.  For  a specific  portion  of  the  structure 
where  additional  benefits  may  be  derived  from  a shape 
redesign,  a substructure  can  be  created  with  the  pertinent 
boundary  load  information  obtained  from  the  FEM.  The  shape 
optimization  problem  for  redesigning  the  shape  of  this 
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substructure  may  then  use  BEM  based  optimization  techniques 
discussed  in  this  study. 

The  current  research  primarily  focused  on  two- 
dimensional  static  structural  shape  optimization  involving 
the  use  of  BEM.  More  recently,  advances  of  the  BEM  in 
problems  related  to  wave  propagation,  structural  vibrations, 
and  three-dimensional  structural  analysis  have  been  reported 
in  [3,44].  The  methodology  developed  in  this  study  for 
efficient  BEM-based  shape  synthesis  can  be  extended  to  these 
disciplines  with  the  attendant  advantages  discussed  in 
earlier  sections. 
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APPENDIX 


In  most  pratical  applications  the  designer  is  interested 
not  only  in  stresses  inside  the  body  and  tractions  on  the 
surface,  but  also  in  the  complete  stress  tensor  on  the 
boundary.  The  calculation  of  this  stress  tensor  proceeds  by 
expressing  the  tractions  in  a local  coordinate  system,  and 
employing  the  strain-displacement  along  a direction,  tangent 
to  the  boundary. 

Assume  a two-dimensional  boundary  element  problem  with  a 
local  Cartesian  coordinate  system  with  axes  tangential  and 
normal  to  the  boundary  (Figure  A.l). 

The  relation  between  the  global  coordinate  x and  local 
coordinate  x is  as  follows. 


The  above  relationship  can  be  expressed  in  terms  of  a 
transformation  matrix  [T]  as  follows. 

x = [T]  x 


For  a plane  strain  problem  in  isotropic  linear 
elasticity,  the  stress  tensor  on  the  boundary  can  be 
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Fingure  A.l  General  two-dimensional  boundary  element 

and  local  coordinate  system  defined  at  the 
stress  point. 
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expressed  in  term  of  tractions  and  the  interior  unit 
extension  tangential  to  the  boundary,  as  follows. 

-11  - Pi 
-12  = P2 

*22  = — — ("  -ii  + 2G"22) 

Hence,  the  stress  tensor  may  be  obtained  as  follows. 

[-]  = [T]T  [a]  [T] 

Using  the  two-dimensional  plane  strain  relations 
and  the  expression, 

[P]  = [T]  [p] 

the  stress  tensor  along  the  boundary  can  be  expressed 
follows . 
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(7  = ( COSa  sin2  a - COS3  a ) Pi 
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